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Abstract 

Structural  biological  materials  such  as  bone,  nacre,  insect  cuticle,  and  sea  sponge 
exoskeleton  showcase  the  use  of  inferior  building  blocks  like  proteins  and  minerals  to 
create  structures  that  afford  load-bearing  and  armor  capabilities.  Many  of  these  are 
composite  structures  that  possess  the  best  of  the  properties  of  their  base  constituents. 
This  is  in  contrast  to  many  engineering  materials,  such  as  metals,  alloys,  ceramics 
and  their  composites  which  show  improvement  in  one  mechanical  property  (e.g.  stiff¬ 
ness)  at  the  cost  of  another  disparate  one  (e.g.  toughness).  These  excellent  design 
examples  from  biology  raise  questions  about  whether  similar  design,  and  improve¬ 
ment  in  disparate  properties,  can  be  achieved  using  common  engineering  materials. 
The  identification  of  broad  design  principles  that  can  be  transferred  from  biological 
materials  to  structural  design,  and  the  analysis  of  the  utility  of  these  principles  have 
been  missing  in  literature.  In  this  thesis,  we  have  firstly  identified  certain  universal 
features  of  design  of  biological  structures  for  mimicking  with  engineering  materials: 
a)  presence  of  geometric  design  at  the  nanoscale,  b)  the  use  of  mechanically  infe¬ 
rior  building  blocks,  and  c)  the  use  of  structural  hierarchies  from  the  nanoscale  to 
the  macroscale.  We  firstly  design,  in  silico,  metal-matrix  nanocomposites,  mimick¬ 
ing  the  geometric  design  found  at  the  nanoscale  in  bone.  We  show  this  leads  to 
improvements  in  flow  strength  of  the  material.  A  key  finding  is  that  limiting  val¬ 
ues  of  certain  of  these  parameters  shuts  down  dislocation- mediated  plasticity  leading 
to  peak  in  flow  strength  of  the  structure.  Metals  are  however,  costly  constituents, 
and  we  next  confront  the  issue  of  whether  it  is  possible  to  use  a  single  mechanically 
inferior  and  commonly  available  constituent,  such  as  silica,  to  create  superior  bio¬ 
inspired  structures.  Wc  turn  to  diatom  cxoskclctons,  protective  armor  structures  for 
algae  made  almost  entirely  of  silica,  and  create  nanoporous  silica  structures  inspired 
from  their  geometry.  We  show  large  improvements  in  ductility  of  silica  through  this 
design,  facilitated  by  a  key  size-dependent  brittle-to-ductile  deformation  transition  in 
these  structures.  Nanostructuring,  while  improving  ductility,  affects  the  stiffness  of 
these  structures,  softening  them  by  up  to  90%  of  bulk  silica.  Hierarchical  assembly 
of  silica  structures  is  then  used  to  regain  the  stiffness  lost  due  to  nanostructuring 


while  not  losing  their  improvement  in  toughness.  Finally,  improvement  in  toughness 
with  several  levels  of  hierarchy  is  studied,  to  showcase  a  defect- tolerant  behavior  that 
arises  with  the  addition  of  hierarchies,  i.e.,  tolerance  of  the  fracture  strength  to  a 
wide  range  of  sizes  of  cracks  present  in  the  structure.  The  importance  of  R-curve 
behavior,  i.e.,  toughness  change  with  the  advance  of  a  crack  in  the  structure,  to  the 
defect-tolcrancc  length  scale  is  also  established.  These  findings  showcase  the  valid¬ 
ity  of  using  design  principles  obtained  from  biological  materials  for  improvement  in 
mechanical  properties  of  engineering  materials. 

Thesis  Supervisor:  Markus  ,1.  Buchler 

Title:  Associate  Professor  of  Civil  and  Environmental  Engineering 


Acknowledgments 


I  would  like  to  take  this  opportunity  to  express  my  heartfelt  gratitude  to  many  who 
have  supported  me  throughout  my  graduate  studies  at  MIT.  First.  I  am  indebted  to 
my  advisor.  Professor  Markus  Buehler,  for  his  mentorship,  advice  and  support  during 
my  PhD.  Discussions  and  brain-storming;  sessions  with  him  have  been  a  great  source 
of  learning  and  inspiration  for  me.  I  am  also  grateful  to  my  peers  in  my  research  group 
for  their  friendship,  discussions,  and  collaborations;  Sinan  Ketcn,  Theodor  Ackbarow, 
Zhiping  Xu.  Andre  Garcia,  Jcrcmie  Bcrtaud,  Steve  Cranford,  Zhao  Qian,  Raffaella 
Paparcone,  Graham  Bratzel,  and  Rouzbeh  Shahsavari.  I  would  like  to  thank  my  col¬ 
laborators  from  various  labs  and  universities:  Prof.  Christian  Thaulow.  Prof.  Pedro 
Reis,  and  Dr.  Kostya  Novoselov.  Brainstorming,  discussing  research  and  writing  pa¬ 
pers  with  them  have  shown  me  the  true  value  of  collaborative  research.  I  would  also 
like  to  thank  my  UROPs  and  other  undergraduate  assistants,  working  with  whom  on 
short  projects  has  been  a  great  teaching  and  learning  experience  for  me. 

The  rigor  of  a  doctoral  work  could  not  have  been  undertaken  without  the  con¬ 
stant  support  of  my  friends  here.  I  am  grateful  to  my  roommates,  Sukant  Mittal, 
Srikanth  Patala  and  Vivck  Inder  Sharma  who  have  stuck  by  me  through  thick  and 
thin.  Coffee  breaks,  and  discussions  with  Srujan  Linga,  Deep  Ghosh,  Vivek  Ragu- 
nathan  and  Sumeet  Kumar  on  every  topic  under  the  sun  have  always  helped  in  breaks 
from  research  work. 

Finally,  I  express  my  deepest  gratitude  to  my  family,  who  have  supported  and 
encouraged  me  in  all  my  academic  endeavors.  I  am  humbled  by  the  sacrifices  made 
by  my  parents,  so  that  I  could  have  a  better  education  and  access  to  world-class 
opportunities. 

This  research  was  funded  by  grants  by  the  Army  Research  Office  (W911NF-06-1- 
0291),  as  well  as  the  UROP  office  at  MIT.  Their  support  is  greatly  appreciated. 

Dedicated  to  my  parents.  Sukla  Sen  and  Siba  Prasad  Sen. 


5 


G 


Contents 


1  Background  19 

1.1  Optimizing  mechanical  properties:  learning  from  nature .  19 

1.2  Structural  hierarchies  m  nature .  22 

1.3  Effect  of  hierarchies  on  failure  across  length  scales .  25 

1.4  Aim  of  study:  hypotheses  .  28 

1.5  Approach  .  29 

1.6  Outline .  29 

2  Methodology  31 

2.1  Atomistic  modeling .  31 

2.1.1  Classical  molecular  dynamics  .  32 

2.1.2  Force  fields  .  35 

2.1.3  Scaling  and  computational  issues .  41 

2.2  Multiscale  modeling .  42 

2.3  Link  to  continuum  state  variables .  44 

2.3.1  Stress  .  44 

2.3.2  Strain .  44 

2.4  Visualization  and  data  analysis  .  45 

2.4.1  Energy  method .  46 

2.4.2  Slip  vector  analysis .  46 

2.4.3  Ccntrosymmetry  parameter .  46 

2.4.4  Common  neighbor  analysis .  47 

2.4.5  Visualization  programs .  48 


7 


3  Strength  enhancement  through  bone-inspired  metal-matrix  nanocom¬ 
posite  design  49 

3.1  Structure  inspiration  from  nanostructure  of  bone .  50 

3.2  Model  construction .  51 

3.3  Interatomic  potential  development  and  testing .  53 

3.3.1  Interatomic  potential  properties .  55 

3.3.2  Design  Parameters .  57 

3.4  Atomistic  simulations  under  tensile  deformation .  58 

3.4.1  Effect  of  geometric  parameters .  58 

3.4.2  Sensitivity  analysis  of  design  parameters .  65 

3.5  Strength  saturation  with  size .  65 

3.6  The  size  confinement  effect  on  dislocation  plasticity .  67 

3.7  Theoretical  analysis  of  the  size  confinement  effect .  68 

3.8  Atomistic  simulations  of  the  size  confinement  effect .  73 

3.9  Results  and  Discussion .  75 

3.10  Conclusions .  83 

4  Ductility  enhancement  through  diatom-inspired  nanoporous  silica 

design  87 

4.1  Background  on  nanoscale  silica  structures  .  88 

4.2  Design  parameters  of  the  nanoporous  silica  structures .  90 

4.3  Materials  and  methods .  90 

4.4  Deformation  of  nano-honeycomb  silica  structures  .  92 

4.5  Analysis  of  deformation  using  theoretical  models .  94 

4.6  Results  and  discussion .  94 

4.6.1  Elasticity  .  96 

4.6.2  Plasticity  and  failure .  104 

4.7  Conclusions .  106 

5  Mesoscale  Model  of  Deformation  and  Failure  of  Hierarchical  Silica 

Nanocomposites  109 


8 


5.1  Review  of  structural  bio-silica  materials .  109 

5.2  Materials  and  Methods .  11 2 

5.2.1  Mesoscale  method  development  and  validation .  112 

5.2.2  Fracture  property  characterization .  116 

5.2.3  R-ciirve  calculation .  117 

5.3  Results  and  discussion .  117 

5.4  Summary  and  Conclusions .  125 

6  The  role  of  multiple  structural  hierarchy  levels  in  defect  tolerance 

and  toughness  129 

6.1  Background  on  structures  with  multiple  levels  of  structural  hierarchy 

and  study  of  their  mechanics .  130 

6.2  Mesoscale  simulation  results .  132 

6.2.1  Two-hierarchy  level  structures  with  periodic  geometry .  132 

6.2.2  Extension  to  3-  and  4-hierarchy  level  structures .  137 

6.3  Discussion  and  Conclusions  .  144 

7  Conclusion  149 

7.1  Summary  of  key  findings  and  significance .  149 

7.2  Opportunities  for  future  research  .  153 


9 


List  of  journal  publications 


1.  D.  Sen,  A.  P.  Garcia,  and  M.  J.  Buehler,  "Mechanics  of  nano-honeycomb  silica 
structures:  A  size-dependent,  brittle-to-ductile  transition",  in  review. 

2.  A.  P.  Garcia,  D.  Sen,  and  M.  J.  Buehler,  “Hierarchical  silica  nanostructures 
inspired  by  diatom  algae  yield  superior  deformability.  toughness  and  strength" , 
Metallurgical  and  Materials  Transactions  A,  accepted. 

3.  D.  Sen,  and  M.  j.  Buehler.  “Atomistically-informed  mesoscale  model  of  defor¬ 
mation  and  failure  of  hierarchical  silica  nanocomposites’5 ,  International  Journal 
of  Applied  Mechanics ,  2010,  2(4). 

4.  D.  Sen,  C.  Thaulow,  S.  V.  Schieffcr.  A.  Cohen,  and  M.  ,1.  Buehler,  “Atom¬ 
istic  Study  of  Crack-Tip  Cleavage  to  Dislocation  Emission  Transition  in  Silicon 
Single  Crystals”,  Physical  Review  Letters,  2010.  104:  p.  235502. 

5.  D.  Sen,  K.  Novoselov,  P.  Reis,  and  M.  J.  Buehler,  “Tearing  graphene  sheets 
from  adhesive  substrates  produces  tapered  nanoribbons”.  Small  2010,  6(10): 
p.  1108.  (Cover  article) 

6.  R.  .lack,  D.  Sen,  and  M.  .1.  Buehler,  “Graphene  nanocutting  through  nanopat.- 
terned  vacancy  defects” ,  Journal  of  Computational  and  Theoretical  Nanosciencc, 
2010,  7:  p.  354-359. 

7.  D.  Sen,  and  M.J.  Buehler,  “Size  and  Geometry  Effects  on  Flow  Stress  in  Bioin¬ 
spired  de  novo  Metal-matrix  Nanocomposites”.  Advanced  Engineering  Materi¬ 
als,  2009,  11(10):  p,  774. 

8.  S.  Cranford,  D.  Sen,  and  M.  J.  Buehler,  “Meso- Origami:  Folding  Multilayer 
Graphene  Sheets".  Applied  Physics  Letters,  2009,  95:  p.  123121. 

9.  I .  Ackbarow,  D.  Sen.  C.  Thaulow  and  M.  j.  Buehler,  “Alpha-Helical  Protein 
Networks  Are  Self-Protective  and  Flaw-Tolerant”,  PLoS  ONE,  2009,  4(0):  p. 
e6015. 


10 


10.  D.  Sen.  and  M.J.  Buehler,  ''Crystal  size  controls  deformation  mechanism: 
Breakdown  of  dislocation  mediated  plasticity  at  uanoscale” ,  Physical  Review 
B,  2008.  77:  p.  195439. 

11.  D.  Sen  and  M.J.  Buehler,  “Shock  loading  of  bone-inspired  metallic  nanocom¬ 
posites".  Solid  State  Phenomena,  2008,  139:  p.  11-12. 

12.  M.  j.  Buehler.  A.  Cohen,  and  D.  Sen,  “Multi-paradigm  modeling  of  fracture 
of  a  silicon  single  crystal  under  mode  II  shear  loading",  .Journal  of  Algorithms 
and  Computational  Technology ,  2008.  2(2):  p.  203-221. 

13.  D.  Sen  and  M.J.  Buehler,  “Simulating  chemistry  in  mechanical  deformation  of 
metals”,  International  Journal  for  Multiscale  Computational  Engineering,  2007. 
5(3-4):  p.  181-202. 


11 


12 


List  of  Figures 


1-1  Toughness-stiffness  Ashby  charts  for  biological  and  engineering  materials.  20 


1-2  Structural  hierarchies  in  bone  and  diatom  exoskeletons .  24 

1-3  Effect  of  hierarchies  on  fracture  resistance  across  scales .  26 

1- 4  Overview  of  computational  methods  across  scales .  30 

2- 1  Molecular  dynamics  simulation  approach .  33 

2- 2  Developing  reactive  potentials .  40 

3- 1  Bone  ultrastructure  and  its  metal-matrix  composite  analogue .  52 

3-2  Properties  of  interatomic  potentials  used .  56 

3-3  Effect  of  volume  fraction  on  flow  stress . 59 

3-4  Effect  of  platelet  offset  on  flow  stress .  60 

3-5  Effect  of  platelet  aspect  ratio  on  flow  stress .  61 

3-6  Effect  of  axial  platelet  spacing  wx  on  flow  stress .  62 

3-7  Effect  of  transverse  platelet  spacing  wy  on  flow  stress .  63 

3-8  Effect  of  interfacial  strength  on  flow  stress .  64 

3-9  Sensitivity  of  flow  stress  to  different  design  parameters .  66 

3-10  Configuration  of  a  thin  strip  geometry  under  shear  loading .  69 

3-11  Morse  and  EAM  potentials  used  for  the  thin  strip  study. .  71 

3-12  Illustration  of  three  regimes  of  slip  events .  72 

3-13  Variation  of  slip  vector  magnitude  on  the  slip  plane .  75 

3-14  Variation  of  ^-component  of  slip  vector  during  the  first  slip  event.  .  .  77 

3-15  Variation  of  ./--component  of  slip  vector  during  the  second  slip  event.  .  78 


13 


3- 16  Application  of  theoretical  model  of  size  effect  to  bone-inspired  nanocom¬ 

posite  behavior .  80 

4- 1  Geometry  for  the  nano-honeycomb  silica  structures  with  different  wall 

widths .  91 

4-2  Stress-strain  graphs  for  silica  nano-honeycomb  structure  for  different 

sizes .  92 

4-3  Von  Miscs  stress  fields  for  silica  nano-honeycomb  structure  for  different 

sizes .  93 

4-4  Geometric  classification  and  stress-strain  behavior  for  additional  nano- 

honeycomb  silica  structures .  95 

4-5  Strain  plots  and  deformation  shape  analysis  for  a  nano-honeycomb 

structure  with  thick  struts .  97 

4-6  Load  distribution  and  stress  transfer  for  the  nano-honeycombs  with 

thick  struts .  98 

4-7  Variation  in  pore  shape  in  the  nano-honeycombs .  100 

4-8  Deformation  shape  analysis  for  a  nano-honeycomb  structure  with  slen¬ 
der  struts .  101 

4-9  Load  distribution  and  stress  transfer  for  the  nano-honcvcombs  with 

slender  struts .  102 

4- 10  Plot  of  yield  stress  versus  a  geometry  parameter  for  nano-honeycomb 

structures  showing  plastic  flow .  105 

5- 1  Structural  hierarchies  in  a  silica-based  skeletal  structures  in  a  sea  sponge.  110 

5-2  Comparison  of  atomistic  and  mesoscale  approaches  for  silica .  113 

5-3  Geometry  of  randomly  distributed  fiber-composite  structures  at  the 

mesoscale .  llg 

5-4  Elastic  moduli  of  the  randomly  distributed  fiber-composite  structures.  119 

5-5  Stress-strain  curves  of  the  randomly  distributed  fiber-composite  struc¬ 
tures .  120 


14 


5-6  Crack  pathways  for  composite  structures  with  nano-honeycomb  struc¬ 
ture  as  the  matrix  and  brittle  silica  as  the  reinforcing  fiber  phase.  .  .  121 

5-7  Crack  pathways  for  composite  structures  with  brittle  silica  as  the  ma¬ 
trix  and  nano-honcycomb  structure  as  the  reinforcing  fiber  phase.  .  .  122 

5-8  Different  composite  structures  with  brittle  bulk-silica  as  the  matrix 

and  ductile  nano-honcycomb  structures  as  the  reinforcing  fiber  phase.  123 
5-9  Calculation  of  the  J-integral  and  R-curves .  124 

5- 10  Effect  of  design  of  the  hierarchical  silica  composite  on  toughness  and 

stiffness  of  bulk  silica .  127 

6- 1  Bone  and  bio-calc-ite  design  morphologies  and  transfer  to  mososcalc 

model .  134 

6-2  Stress-strain  curves  of  the  bone-like  and  biocalcite-like  nanoporous  sil¬ 
ica/  bulk  silica  structures .  135 

6-3  Crack  pathways  for  composite  structures  with  bone-like  and  biocalcite- 

like  desi  gn .  136 

6-4  Self-similar  (fractal)  and  dissimilar  3- hierarchy  structures .  138 

6-5  Comparison  of  stress-strain  behavior  of  the  2-hicrarchy  and  3-hicrarehy 

self-similar  nanoporous  silica/  bulk  silica  structures .  139 

6-6  Comparison  of  stress-strain  behavior  of  the  2-hierarchy  and  3-hierarchv 

dissimilar  nanoporous  silica/  bulk  silica  structures .  140 

6-7  Mechanisms  of  defect-tolerance  in  the  3-hierarchy  dissimilar  structures.  142 
6-8  4-hierarchy  structure  morphology,  their  stress-strain  plots  and  R-curve 

behavior .  143 

6-9  Link  between  R-curve  behavior  and  unstable  crack  propagation.  .  .  .  145 

6-10  Loss  of  fracture  strength  vs.  crack  size .  146 

6-11  Defect- tolerant  length  scale  vs.  the  number  of  hierarchies .  147 


15 


16 


List  of  Tables 


3.1  Work  of  adhesion  across  a  “strong"  and  "weak”  interface .  55 

3.2  Summary  of  the  critical  length  scales  for  shutdown  of  dislocation  plas¬ 
ticity  in  single  crystal  Cu  and  A1 .  72 

4.1  Elastic  moduli  for  liano-lioneycomb  structures  with  thick  struts.  .  .  .  99 

4.2  Elastic  moduli  for  nano-honeycomb  structures  with  slender  struts.  .  .  103 


17 


18 


Chapter  1 


Background 

1.1  Optimizing  mechanical  properties:  learning  from 
nature 

There  exist  several  natural  materials  with  exceptional  mechanical  properties.  Silk  has 
a  strength  per  unit  weight  larger  than  steel;  bone,  nacre  and  sea  shell  have  excellent 
toughness  properties  given  their  weak  and  brittle  constituents  (hydroxyapatite,  silica, 
protein ) .  Several  of  these  structural  biological  materials  are  composites,  with  very 
good  deformation-resisting  and  load-carrying  capacities.  The  composite  constituents 
are  usually  different  proteins,  such  as  collagen  or  chitin.  and  minerals  such  as  calcitc, 
aragonite  and  hydroxyapatite.  These  composites  arc  typically  lightweight,  but  possess 
an  unusual  stiffness,  strength,  toughness  and  fatigue  resistance  for  their  composition 
and  weight.  Very  importantly,  the  final  material  properties  are  seen  to  be  a  combi¬ 
nation  of  the  best  properties  of  the  base  constituents  and  cannot  be  approximated  a 
rule-of-mixtures  calculations  of  their  base  material  properties  [1].  Most  engineering 
analogs,  on  the  other  hand,  show  the  characteristic  ‘banana-curve’  type  behavior, 
that  is,  they  do  not  enable  the  combination  of  high  levels  of  strength,  stiffness  and 
toughness  (Figure  1-1). 

There  is  nothing  remarkable  about  the  mechanical  properties  of  the  individual 
constituents  of  these  natural  materials.  Hydroxyapatite  or  silica  possess  fracture 
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Figure  1-1:  (a)  Toughness  versus  stiffness  for  a  number  of  biological  materials  (based  on 
the  data  compilation  in  [2].  Biological  composites,  such  as  antler,  dentin,  bone  and  enamel 
combine  the  good  properties  of  the  protein  and  mineral  components  and  are  typically  both 
stiff  and  tough.  Generally  speaking,  the  stiffness  increases  and  the  toughness  decreases  with 
mineral  content  from  antler  to  dentin/bone  and  enamel  [3];  (b)  Toughness  vs.  stiffness  for 
metals,  alloys,  ceramics,  and  metal-matrix  composites  (MMC),  lie  on  the  ‘banana-curve’,  an 
inverse  relation  between  increasing  toughness  and  decreasing  stiffness.  The  circular  yellow 
region  shows  the  property  region  of  high  toughness  and  stiffness  that  may  be  accessible 
through  designing  bio-inspired  composite  structures.  Figure  adapted  from  [4], 
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toughness  much  like  man-made  ceramics,  proteins  have  stiffness  close  to  synthetic 
polymers.  It  is  thus  to  be  reasoned  that  it  is  the  structure  and  geometry  of  these 
materials  that  give  rise  to  excellent  properties.  A  key  common  feature  of  these  com¬ 
posite  biological  materials  seen  is  the  presence  of  remarkable  designs  with  building 
blocks,  often  hierarchically  arranged  from  the  nanometer  to  the  macroscopic  length 
scales  [5.  6.  7]  that  give  rise  to  excellent  properties.  Every  structural  level  (hierarchy) 
is  postulated  to  contribute  to  the  mechanical  stability  and  properties  of  the  resulting 
design.  This  endows  them  with  a  far  greater  level  of  structural  complexity  and  or¬ 
ganization  than  synthetic  composites.  However,  the  extent  to  which  the  presence  of 
hierarchies  quantitatively  affects  structural  properties  is  not  clearly  understood,  and 
a  broad  framework  for  studying  structural  hierarchies  is  missing  in  prior  literature. 
Hierarchical  organizations  have  already  been  studied  in  numerous  fields  far  removed 
from  structural  materials,  c.g.  ecology,  transportation,  and  engineering  control  sys¬ 
tems  [8,  9,  10],  and  recent  attempts  to  draw  parallels  from  these  fields  to  the  study 
of  hierarchies  in  biological  structures  point  to  exciting,  uncharted  territories  in  the 
study  of  structure  of  biological  materials  [11,  12,  13]. 

An  application  of  this  emerging  science  would  be  in  the  design  of  synthetic  hierar¬ 
chical  materials  [14].  The  materials  used  can  be  based  on  metals,  ceramics,  polymers 
and  proteins  [15,  16,  17,  18].  The  broad  vision  would  be  to  design  these  materials  us¬ 
ing  a  bottom-up  approach,  starting  at  the  nanoscalc,  and  building  in  structure,  as  the 
length  scale  is  increased.  The  design  at  the  lowest  length  scales  (nano-  to  micro-)  is 
particularly  suited  for  computational  nanoscale  experiments,  as  carried  out  through 
molecular  dynamics  (MD)  and  mesoscale  simulations.  Atomic  scale  molecular  dy¬ 
namics  simulations  have  been  used  previously  to  probe  size-scale  and  temperature 
scale  effects  in  the  deformation  of  nanomaterials,  often  providing  design  templates 
[19.  20,  21]. 

Many  of  these  biological  materials  are  seen  to  possess  the  following  design  features: 
a)  use  of  inferior  base  constituent  materials,  (b)  composite  sub-structure  starting  at 
nanoscale  dimensions;  and  (c)  hierarchical  arrangement  of  structure  from  nanoscale 
to  macroscale.  We  describe  what  hierarchies  are  in  the  next  section  and  provide  a 
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few  examples. 


1.2  Structural  hierarchies  in  nature 

Structural  hierarchies  in  natural  biological  materials  are  defined  as  the  quality  of 
certain  materials  possessing  structure  and  organization  at  several  levels  of  length 
scale  [22.  23].  A  natural  definition  of  hierarchical  systems  arises  from  the  field  of 
systems  theory  [24.  25],  where  they  are  defined  as  composition  of  stable,  observ¬ 
able  sub-elements  that  are  unified  by  a  super-ordinate  relation  [8].  The  stability  of 
sub-elements  at  different  levels  makes  them  building  blocks  for  the  next  higher  level. 
Averaging  or  coarse-graining  of  properties  over  one  hierarchy  level  to  derive  informa¬ 
tion  for  the  next  higher  level  is  usually  not  feasible,  owing  to  linking  of  behavior  across 
several  levels,  unless  there  is  a  large  separation  of  length  scales  across  successive  levels 
[22]. 

A  key  feature  of  hierarchical  biological  systems,  structural  or  otherwise,  is  their 
robustness.  Robustness  has  been  studied  by  systems  scientists  for  biological  systems 
and  classified  in  the  following  ways:  (a)  adaptation  -  the  ability  to  cope  with  external 
changes,  (b)  parameter  insensitivity,  and  (c)  graceful  degradation-  slow  degradation 
of  a  system's  function  after  damage,  rather  than  catastrophic  failure  [25].  Hierarchical 
systems  in  biology,  studied  from  a  system-theoretic  point  of  view  show  optimality  of 
several  properties  and  robustness,  at  the  same  time  [11],  Hierarchical  systems  also 
show  improved  behavior  (optimality)  over  a  large  number  of  mechanical  properties 
simultaneously,  as  compared  to  their  base  elements  at  the  lowest  level  of  hierarchy. 
Whether  individual  properties  are  improved  by  design  elements  at  individual  levels 
of  hierarchy,  or  arise  from  a  combination  of  properties  at  different  scales  is  debatable 
(e.g.  sec  Section  1.3,  “Effect  of  hierarchies  on  failure  across  length  scales’’). 

Here,  we  provide  two  examples  of  hierarchical  structural  geometries.  Human  corti¬ 
cal  bone  [26]  is  seen  to  be  composed  of  7  structural  levels  of  hierarchical  arrangement 
and  possesses  excellent  strength  and  toughness  properties  while  being  lightweight. 
This  is  thus  an  ideal  material  to  mimic  for  high  toughness  applications,  while  probing 
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effect  of  hierarchies  on  mechanical  properties.  Bone  is  a  composite  of  organic  and 
inorganic  constituents:  30%  bone,  by  weight,  is  organic:  of  which  90-95%  is  collagen, 
rest  is  non-collagenous  proteins.  At  the  nanoscale,  bundles  of  collagen  molecules  are 
arranged  in  fibrils,  which  are  twisted  in  a  coil  (fiber).  70%  of  bone  is  made  up  of  the 
inorganic  mineral  hydroxyapatite,  which  includes  calcium  phosphate,  calcium  car¬ 
bonate,  calcium  fluoride,  calcium  hydroxide  and  citrate.  This  inorganic  component 
(([Ca3(PO)4)2]3  •  Ca(OH)2))  is  predominantly  crystalline.  The  crystals  are  platelets 
or  rods,  about  8  to  15  A  thick,  20  to  40  A  wide  and  200  to  400  A  long  and  arranged 
in  a  regular  array  at  the  nanoscale  (see  Figure  1-2).  The  mineralized  collagen  fibers 
form  planar  arrangements  called  lamellae  (3-7  pm  wide).  These  sheets  (lamellae)  of 
mineralized  collagen  fibers  wrap  in  concentric  layers  around  a  central  canal  to  form 
osteons.  Osteons  appear  like  cylinders  ~200-250  pm  in  diameter  running  parallel  to 
the  long  axis  of  the  bone.  Figure  1-2  shows  all  these  levels  of  hierarchy  in  bone  from 
the  nanoscale  up  to  the  macroscale. 

Diatom  exoskeletons  are  a  composite  of  97%  silica  and  remaining  3%  protein  mate¬ 
rial  such  as  silieat.eins.  The  structure  consists  of  4  levels  of  hierarchy.  At  the  nanoscale 
is  the  basic  constituent,  biosilica,  made  of  fused  silica  nanospheres  connected  via  an 
organic  matrix.  This  biosilica  is  designed  into  a  porous  nanostructure,  the  cribellum, 
at  the  smallest  assembly  scale,  in  a  regular  lattice  arrangement  of  pores  with  pore 
diameters  ~45  nm  and  distances  ~68  nm.  The  second  layer,  the  cribrum  possesses 
a  hexagonal  porous  structure  with  larger  pore  sizes  of  «200  nm.  The  largest  porous 
layer,  the  areola  has  a  pore  size  of  «1.1  pm  [27]. 

The  importance  of  hierarchical  design  to  improving  divergent  mechanical  prop¬ 
erties  in  biological  applications  has  been  proposed  by  many  authors.  For  example, 
the  structural  hierarchy  seen  in  skeleton  of  sea  sponge  [28,  29,  30]  is  supposed  to 
be  responsible  for  its  high  strength  and  crack  resistance  despite  being  made  almost 
completely  of  brittle  silica.  Hierarchical  arrangements  in  protein  structures  from 
amino  acids  up  to  secondary  structures,  have  been  proposed  as  an  arrangement,  for 
improving  robustness  of  the  structures  [13].  Hierarchical  assembly  thus  might  hold 
the  key  to  scaling  up  excellent  mechanical  properties  seen  in  many  synthetic  nanonia- 
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Figure  1-2:  Structural  hierarchies  in  two  different  biological  materials,  (a)  bone,  showing 
from  left  to  right,  cortical  and  cancellous  bone  (different  types  of  bone);  osteons;  lamellae; 
collagen  fiber  assemblies  of  collagen  fibrils;  bone  mineral  crystals,  collagen  molecules,  and 
non-collagenous  proteins  (Figure  reproduced  from  [1]);  and  (b),  marine  diatom  species 
( Concinodicus  sp.),  a  silica-based  cxoskclcton  (called  frustule)  made  up  of  porous  parts 
arranged  in  a  hierarchical  fashion,  showing,  from  left  to  right,  the  whole  frustule  member 
(external  surface  of  the  diatom);  areola  pores,  the  internal  surface  of  the  diatom;  the  2nd 
central  porous  layer,  the  cribrurn;  the  cribellum,  the  external  porous  layer.  The  three  layers 
are  arranged  on  top  of  each  other  (Figure  reproduced  from  [27]). 
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terials  up  to  macroscale  engineering  structures.  Carbon  nanotubes,  graphene,  metal 
nanowires  have  excellent  strength  but  cannot  be  presently  used  in  engineering  struc¬ 
tures  because  the  procedure  for  connecting  disparate  length  scales  while  maintaining 
nanoscale  properties  remains  unknown  [31], 


1.3  Effect  of  hierarchies  on  failure  across  length 
scales 

Loading  and  fracture  experiments  on  structural  biological  materials  have  revealed 
a  complex  set  of  mechanisms  across  wide  length  scales.  The  question  of  which  are 
the  dominant  mechanisms  in  the  failure  of  a  particular  material,  and,  are  hierarchies 
and  the  multiple  scale  mechanisms  they  engender,  essential  for  the  improved  frac¬ 
ture  properties,  arc  still  hotly  debated.  Experimentally,  this  is  difficult  to  observe 
because  of  the  difficulty  of  testing  substructures  of  a  material  at  different  length 
scales,  and  the  problem  with  separating  the  contribution  of  different  mechanisms  to 
the  overall  toughness.  Here,  we  briefly  review  the  experimental  evidence  of  effect 
of  hierarchies  on  failure  in  two  biological  systems  with  different  basal  components, 
bone  (hydroxyapatite-based)  and  diatoms  (silica-based).  These  materials  are  chosen 
as  two  representative  systems  which  we  will  use  as  design  templates  using  engineering 
materials  in  later  sections  of  the  thesis. 

A  major  property  of  bone  is  its  fracture  resistance  and  toughness.  This  has  been 
attributed  to  distinct  mechanisms  on  different  length  scales  by  various  authors.  On 
the  micron  length  scale,  where  the  bone  structure  consists  of  osteons,  this  lias  been 
attributed  to  two  mechanisms  (a)  crack  bridging  and  (b)  microcracking  [32],  On 
the  nanometer  length  scale,  this  has  been  attributed  to  flaw  tolerance  of  size  of 
mineral  platelets  [33].  The  modular  domain  nature  of  the  organic  matrix  at  different 
scales,  and  its  stepwise  unfolding  has  been  proposed  as  a  mechanism  for  the  intrinsic 
toughness  of  the  protein  matrix  [34].  Below  we  provide  brief  descriptions  of  both 
mechanisms. 
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Figure  1-3:  Toughening  mechanisms  exist  at  different  levels  of  hierarchy  in  bone.  These 
arc  molecular  uncoiling  and  sliding  at  the  protein  domain  level,  and  flaw-tolcrance  at  the 
individual  hydroxyapatite  crystal  level,  microcracking,  crack-bridging,  and  crack  deflection 
at  the  larger  scales.  Figure  reproduced  from  [35]. 
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From  critical  stress  intensity  factor  ( KIc )  studies  of  bending  loading  of  notched 
bone  specimen  to  failure,  formation  of  microcracks  are  seen,  which  are  of  the  order 
of  osteons  in  length.  It  is  hypothesized  that  microcracks  tend  to  originate  around 
osteons  due  to  debonding  at  osteon-matrix  interface  or  osteon  pull-out.  The  presence 
of  microcracks  in  the  wake  of  a  crack  have  been  shown  to  result  in  the  residual  opening 
of  the  crack  tip,  and  a  redistribution  of  stresses  in  the  crack  tip  region,  which  reduces 
the  crack  extension  force  and  increases  the  toughness  of  the  material.  Crack  bridging 
in  the  wake  of  a  crack  has  also  been  proposed  as  a  crack  tip  shielding  mechanism. 
Crack  bridging  involves  formation  of  unbroken  regions  that  span  the  crack  in  the 
wake  of  the  crack  tip  and  act  to  resist  crack  opening.  Such  bridging  can  results  from 
uncracked  ligaments  and  intact  collagen  fibrils.  Both  these  mechanisms  thus  reduce 
crack  propagation,  and  the  dominant  mechanism  is  still  under  debate  [32]. 

Another  research  direction  has  been  considering  the  characteristic  nanostructure 
of  bone,  a  geometric  motif  that  is  common  to  other  hard  structural  biomaterials  such 
as  nacre  and  dentin.  The  nanostructure  of  bone  is  seen  in  Figure  1-2.  consisting  of 
mineral  platelets  arranged  in  a  staggered  pattern  in  a  collagen  matrix.  The  com¬ 
monality  of  this  structural  motif  across  structural  materials  suggests  some  intrinsic- 
properties  in  the  design  that  improve  mechanical  properties.  A  mechanism  has  been 
proposed  by  Gao.  Fratzl  et  al.  [36]  whereby  under  tensile  loading,  staggered  mineral 
platelets  carry  tensile  load  and  the  protein  matrix  transfers  the  load  between  mineral 
crystals  via  shear.  The  fracture  toughness  of  the  composite  depends  on  the  tensile 
strength  of  the  mineral  platelets.  It  has  been  showed  that  the  nanoscale  width  of 
the  mineral  platelets  embedded  in  the  collagen  matrix  is  such  that  the  material  be¬ 
comes  insensitive  to  crack  like  flaws  at  this  length  scale  (approximately  30  inn)  and 
fails  under  tension  at  the  theoretical  strength  for  a  perfect  crystal  [33].  When  the 
mineral  size  exceeds  a  length  scale  of  order  of  30  nm,  fracture  strength  is  sensitive  to 
structural  size.  This  concept  is  called  the  flaw  tolerance  of  size  of  mineral  platelets. 
The  theory  claims  that  the  size  of  mineral  platelets  in  bone  is  optimized  at  this 
flaw  tolerant  size.  The  optimum  aspect  ratio  (height/width)  of  these  platelets  can 
be  obtained  by  assuming  that  protein  and  mineral  fail  at  the  same  time.  However, 
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major  shortcomings  of  this  simplified  model  are  that  it  fails  to  take  into  account  the 
complex  non-stoichiometric  chemistry  at  mineral-protein  interfaces  [37.  38]  and  size 
limitations  of  mineral  owing  to  the  same. 

Diatom  algae  form  a  very  different  kind  of  protective  exoskeleton.  This  protective 
covering  is  porous  and  made  of  up  to  97%  silica.  Several  studies  reported  in  the 
recent  literature  have  revealed  the  mechanical  properties  of  diatom  shells.  Hamm  et 
al.  [39]  used  a  glass  needle  to  load  and  break  diatom  frustules  in  order  to  probe  their 
mechanical  response  at  failure,  and  found  high  strength  (between  «1  and  7  MPa 
compressive  stress  for  fracture)  and  reversible  elastic  strains  (e.g.  2.5%  reversible 
strain  in  a  frustulc  section).  A  three-dimensional  finite-element  model  of  the  frustule 
in  the  same  work  showed  that  the  highest  stresses  within  the  frustule  before  failure  was 
«540  MPa.  The  geometric  design  of  the  frustule  led  to  the  applied  external  pressures 
creating  homogeneous  stress  distributions  within  the  Structure.  Other  researchers 
[40.  41]  have  used  AFM  nanoindentation  to  study  the  nanoscale  material  properties 
of  the  porous  frustule  layers  of  diatoms,  identifying  pore  sizes  on  the  order  of  several 
tens  of  nanometers  at  the  smallest  levels  in  the  hierarchy,  with  ultra-thin  silica  walls 
on  the  order  of  several  nanometers.  They  observed  that  the  variation  of  mechanical 
properties  between  the  hierarchical  frustule  layers  could  be  influenced  by  the  pore 
size,  pore  distance,  porosity,  and  under  different  biomineralization  processes. 

This  brief  review  of  the  physics  of  toughness  and  fracture  strength  of  bone  and 
diatoms  show  the  importance  of  hierarchical  levels  in  optimizing  mechanical  proper¬ 
ties.  Structures  at  very  different  length-scales  have  been  seen  to  play  a  prominent 
role  in  their  stiffness,  strength  and  toughness  [32,  33,  41]. 


1.4  Aim  of  study:  hypotheses 

The  aim  of  this  thesis  is  the  investigation  of  the  effect  of  mimicking  certain  universal 
features  found  in  biological  structural  materials,  and  to  assess  their  potential  for  use 
in  the  design  of  engineering  materials.  The  key  hypotheses  are  that  these  universal 
features  that  can  be  transferred  to  the  design  of  engineering  materials  are: 
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(a)  geometric  design  at  the  ultimate  scale  (nanostructure)  of  biological  structure, 
with  the  underlying  effects  on  deformation  mechanisms  at  the  nanoscale. 


(b)  the  use  of  mechanically  inferior  constituent  materials, 

(c)  attempting  bottom-up  design  with  hierarchical  molecular-scale  (nm)  assemblies 
up  to  macroscale  (mm/cm)  dimensions. 

The  particular  individual  materials  used  in  the  design  is  not  the  important  con¬ 
cept,  rather  it  is  differentiating  the  requirement  of  nanostructuring  and  hierarchical 
assembly  in  the  improvement  of  mechanical  properties.  Addressing  these  issues  by 
laying  out  a  computational  and  theoretical  framework  is  the  fundamental  goal  of  this 
thesis. 


1.5  Approach 

The  scope  of  our  study  lends  itself  to  the  use  of  bottom-up  computational  simulations 
as  a  design  and  analysis  tool.  Figure  1-4  shows  computational  methods  across  size 
and  time  scales  that  arc;  available  in  the  computational  mechanics  literature. 

In  this  thesis,  we  study  the  mechanical  properties  at  the  nanoscale  using  atomistic 
simulations,  and  at.  the  sub-micron  and  micro-scale  using  mesoscale  modeling.  These 
methods  are  outlined  in  greater  detail  in  Chapter  2. 

1.6  Outline 

The  content  of  this  thesis  is  arranged  as  follows:  Chapter  2  outlines  all  of  the  com¬ 
putational  methods  used  in  this  thesis,  from  flic  atomistic  to  nmltisealc  modeling 
approaches.  It  also  provides  some  background  on  the  analysis  techniques,  and  visu¬ 
alization  methods  used.  In  chapter  3,  we  take  inspiration  from  the  design  of  bone 
nanostructure  and  describe  the  design  and  mechanical  properties  of  a  metal-matrix 
nanocomposite  based  on  it.  Metals,  however,  are  costly  materials  to  use  as  con¬ 
stituents.  Is  it,  possible  to  use  cheaper  and  more  readily  available  materials  as  build- 
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Figure  1-4:  Overview  of  computational  methods  across  scales,  scanning  from  quantum 
calculations  at  the  Angstrom  scale  to  continuum  models  at  the  meter  scale.  Multiscale 
coupling  of  methods  can  be  used  to  traverse  through  a  wide  range  of  time-  and  length 
scales.  Figure  adapted  from  [42], 


ing  blocks,  such  as  silica,  abundantly  found  in  sand?  Bulk  silica,  however,  is  a  ‘weak’ 
material  for  structural  purposes,  due  to  low  toughness  and  brittle  failure  behavior. 
We  tackle  the  problem  of  enhancing  the  ductility  of  silica  in  chapter  4,  by  taking 
inspiration  from  a  silica-based  biomaterial,  the  diatom,  and  describe  the  design  and 
mechanical  properties  of  a  nanoporous  silica  structure  inspired  from  its  nanostruc¬ 
ture.  The  next  question  is  how  to  measure  the  mechanical  properties  of  hierarchical 
structures  built  from  the  bottom-up  using  these  nanostructures,  and  we  proceed,  in 
Chapter  5.  to  develop  a  mesoscale  modeling  method  to  describe  the  mechanics  of  hier¬ 
archical  silica  nanocomposites.  Chapter  G  describes  the  use  of  this  mesoscale  method 
in  measuring  toughness  improvements  over  several  levels  of  hierarchy.  Chapter  7 
summarizes  major  findings  of  this  thesis,  and  provides  an  outlook  for  future  research. 
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Chapter  2 


Methodology 


In  this  Chapter,  a  brief  overview  of  computational  techniques  used  in  this  work  is 
presented.  The  focus  is  on  computational  methods  to  study  mechanical  behavior  of 
materials,  and  they  are  classified  here  by  the  length  scales  accessible  to  the  different 
methods.  Firstly,  different  schemes  of  atomistic  simulations  are  reviewed  which  can 
access  nanometer  and  sub-nanometer  length  scales.  Then  we  cover  multiscale  mod¬ 
eling  methods  which  can  access  sub-micron  and  micron  length  scales.  In  the  next 
section,  we  establish  the  link  between  the  data  from  atomistic  and  mcsoscale  simu¬ 
lations  and  material  properties  measured  at  a  continuum  level.  Finally,  visualization 
techniques  and  packages  are  covered. 


2.1  Atomistic  modeling 

This  section  describes  the  atomistic  modeling  approaches  used  in  this  thesis.  Par¬ 
ticular  emphasis  is  given  to  the  molecular  dynamics  method  and  its  variants.  An 
extensive  discussion  of  atomistic  force  fields,  which  are  a  key  factor  in  the  success 
of  any  atomistic  simulation,  is  also  undertaken.  Finally,  scaling  and  computational 
limitations  of  atomistic  modeling  are  briefly  touched  upon. 
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2.1.1  Classical  molecular  dynamics 


Molecular  dynamics  (Ml))  is  a  computer  simulation  tool  for  studying  the  real-time 
motion  of  a  group  of  atoms  or  molecules  under  their  mutual  interactions  for  a  cert  ain 
period  of  time.  Since  its  development  in  the  late  1950s  [43,  44],  it  has  been  used  to 
simulate  groups  of  atoms  in  sizes  from  a  few  to  several  billions  recently,  and  over  a  time 
of  femtoseconds  (fs)  to  fractions  of  microseconds.  It  has  found  several  applications 
in  materials  science,  chemistry,  solid-state  physics,  fluid  mechanics,  biomechanics 
and  other  fields.  In  the  field  of  materials  science,  the  MD  method  is  capable  of 
capturing  atomistic  mechanisms  that  play  a  key  role  in  several  materials  phenomena 
e.g.  deformation,  fracture,  diffusion,  chemical  reactions,  self-assembly,  and  phase 
transformations. 

The  core  requirement  for  an  atomistic  simulation  using  MD  is  a  (2-body  to  multi- 
body)  atom  interaction  potential  or  force  description.  This  is  a  coarser  system  descrip¬ 
tion  than  the  use  of  ab-initio  quantum  methods  [45]  to  describe  atomic  interactions 
that  allows  an  accurate  description  of  ground-state  and  excited-electronic  states  of  a 
group  of  atoms,  leading  to  ab-initio  molecular  dynamics.  However,  ab-initio  methods 
are  hugely  expensive  in  terms  of  computational  time  and  power,  and  can  only  be  used 
for  systems  with  a  few  100  atoms.  What  we  will  discuss  further  is  classical  molecular 
dynamics,  where  the  ground-state  of  the  system  is  used  to  obtain  an  atomic  interac¬ 
tion  potential,  losing  all  electronic  states  information.  The  atomistic  potentials  can 
however  be  designed  to  capture  reactions  and  bond  making  or  breaking  events.  The 
main  idea  then  is  to  compute  the  dynamical  trajectory  of  each  atom  in  the  group, 
considering  their  atomic  interaction  potentials,  by  solving  each  atom's  equation  of 
motion  according  to  F  =  ma,  where  F,  m  and  a  are  force,  mass  and  acceleration 
respectively,  leading  to  atomic  positions  r,(f),  velocities  Vi(t)  and  accelerations  at(t) 
for  the  zth  atom.  The  numerical  integration  of  Newton’s  law  by  considering  proper 
interatomic  potentials  to  obtain  interatomic  forces  enables  one  to  simulate  a  small/ 
large  group  of  atoms.  The  basic  concept  of  molecular  dynamics  is  shown  in  Figure 
2-1  (a, b). 
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Figure  2-1:  Model  of  the  individual  energy  contributions  due  to  bond  stretching,  bond 
bending,  bond  rotation  as  well  as  electrostatic  and  vdW  interactions.  The  combination 
of  these  terms  constitutes  the  entire  energy  landscape  of  interatomic  and  intermolecular 
interactions.  Figure  reprinted  from  Ref.  [46]. 


Classical  molecular  dynamics  generates  the  trajectories  of  a  large  number  of  par¬ 
ticles,  interacting  with  a  specific  interatomic  potential.  Thereby,  the  complex  3D 
structure  of  an  atom  (composed  of  electrons  and  a  core  of  neutrons  and  protons)  is 
approximated  by  a  point  particle,  as  shown  in  Figure  2-  1(b)).  The  total  energy  of 
the  system  is  written  as  the  sum  of  kinetic  energy  (A')  and  potential  energy  (£/), 


E  =  K  +  U 


(2.1) 


where  the  kinetic  energy  is 


K  — 


iV 

E 

3= 1 


9 


and  the  potential  energy  is  a  function  of  the  atomic  coordinates  rr 


(2.2) 


U 


U(rj), 


(2.3) 
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with  a  properly  defined  potential  energy  surface  |/(?q).  The  numerical  problem  to  be 
solved  is  a  system  of  coupled  second  order  nonlinear  differential  equations: 


3  i--v-  (2.4) 

which  can  only  be  solved  numerically  for  more  than  two  particles.  N  >  2.  Typically, 
MD  is  based  on  updating  schemes  that  yield  new  positions  from  the  old  positions, 
velocities  and  the  current  accelerations  of  particles.  In  the  commonly  used  Verlet 
scheme,  this  can  be  mathematically  formulated  as 

r;:(/u  +  At)  =  — r?;(f0  —  At)  +  2rl(t0)  +  aj(to)  (At)-  +  0((At)4)  (2.5) 

The  forces  and  accelerations  are  related  by  a{  =  fi/m.  The  forces  are  obtained 
from  the  potential  energy  surface  -  sometimes  also  called  force  field  -  as 

F  =  =  -  VrJ'ii-j)  j  =  1..N.  (2.6) 

Molecular  dynamics  simulations  can  be  used  as  a  tool  for  discovering  mechanisms 
or  reaction  pathways  in  a  small  group  of  atoms,  where  the  simulated  system  is  the 
actual  designed  experimental  sample.  Thus  it  can  be  used  for  small  molecule  reac¬ 
tions.  and  nano-sized  systems  such  as  carbon  nanotubes,  proteins,  nanofluids  etc. 
However,  it.  is  also  extensively  used  to  create  a  sample  picture  of  much  larger  systems 
in  length  scale.  The  idea  then  is  to  use  molecular  dynamics  as  a  tool  for  statistical 
mechanics  for  studying  systems  under  equilibrium  or  evolving  boundary  conditions. 
In  either  case,  the  actual  system,  may  be  under  certain  thermodynamic  conditions 
e.g.  a  certain  temperature  or  pressure  or  energy  constraint.  These  thermodynamic 
boundary  conditions  need  to  be  applied  to  the  MD  simulation.  Time- averaged  ther¬ 
modynamic  variables  can  then  be  estimated  over  the  length  of  the  MD  run.  The 
ergodie  hypothesis  then  postulates  that  the  time  averaged  statistical  quantities  for  a 
system  approach  the  ensemble  average  over  all  possible  states  of  the  system  at  very 
large  times. 

The  primarily  used  thermodynamics  conditions,  or  statistical  ensembles  in  MD 
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are  the  micro-canonical  (also  called  NVE  -conserved  number  of  atoms,  volume  and 
energy),  the  canonical  (also  called  NVT-  conserved  number  of  atoms,  volume  and 
temperature)  and  the  isobaric-isothermal  (also  called  NPT  -  conserved  number  of 
atoms,  pressure  and  temperature).  Next  we  provide  a  brief  overview  of  how  these 
ensembles  are  implemented  in  various  molecular  dynamics  codes. 

2. 1.1.1  NVE  ensemble 

In  the  NVE  ensemble,  the  group  of  atoms  is  isolated  from  any  changes  in  number  of 
atoms  ( N ),  system  volume  ( V),  and  total  system  energy  (E).  The  system  evolution 
is  by  the  same  equations  provided  in  (2.5).  It  is  critical  to  ensure  energy  conservation 
in  the  numerical  approximation  of  these  equations,  which  require  a  judicious  choice 
of  the  timestep  of  integration  in  (2.5). 

2. 1.1. 2  NVT  ensemble 

In  the  NVT  ensemble,  the  number  of  atoms  (iV),  system  volume  ( V)  and  system 
temperature  (T)  are  conserved.  The  system  is  thus  allowed  to  exchange  energy 
with  a  virtual  heat  bath,  to  maintain  constant  temperature.  The  implementation  is 
through  the  interactions  of  the  atoms  in  the  system  with  a  thermostat.  A  simple  and 
common  implementation  is  the  Berendsen  thermostat  [47],  which  scales  the  velocities 
of  the  atoms  every  few  steps  in  the  simulation  so  that  the  temperature  approaches  the 
desired  value,  thus  mimicking  a  heat  bath.  This  is  realized  by  calculating  a  rescaling 
parameter  A , 


Wi+f  (£-i)-  (2-7) 

where  At.  is  the  MD  time  step  and  r  is  a  parameter  that  describes  the  strength  of  the 
coupling  of  the  system  to  the  virtual  heat  bath.  The  velocities  are  then  rescaled  by, 

Viea'.j  (2.8) 
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for  each  atom  i,.  Other  approaches  to  enforce  the  NVT  ensemble  include  the  Nose- 
Hoover  scheme  [48]  and  methods  based  on  Langevin  dynamics  [49]. 


2. 1.1. 3  NPT  ensemble 

In  the  NPT  ensemble,  the  number  of  atoms  (N),  system  pressure  (P)  and  system 
temperature  (2’)  are  conserved.  In  addition  to  a  thermostat,  a  barostat  is  needed, 
and  the  system  volume  is  adjusted  for  the  system  pressure  to  converge  to  an  applied 
pressure  tensor.  Hein  the  popular  schemes  are  the  Nose-Hoover  [50],  and  Parrinello- 
Rahman  [51] . 

The  availability  of  interatomic  potentials  for  a  specific  material  is  often  a  limiting 
factor  for  the  applicability  of  the  MD  method,  since,  the  complete  material  behavior  is 
determined  by  this  choice.  Designing  appropriate  models  for  interatomic  interactions 
provides  a  rather  challenging  and  crucial  step  that  remains  the  subject  of  very  active 
discussions  in  the  scientific  community.  A  variety  of  different  interatomic  potentials 
are  used  in  the  studies  of  metals,  inorganic  materials  and  biological  materials,  with 
various  degrees  of  accuracy  and  speed,  and  the  choice  of  potential  depends  heavily  on 
the  application  area,  and  not  just,  a  stand-alone  description  of  the  material  behavior. 
Designing  or  choosing  a  potential  that  is  accurate  enough  to  capture  all  possible 
phenomena  expected  in  the  application  under  study,  and  not  too  complex  that  it 
slows  down  computation  time  and  size  scale  considerably  is  an  art  in  this  field. 

The  goal  of  the  next  section,  2.1.2.  is  to  provide  a  brief  overview  of  popular 
interatomic  force  fields  and  modeling  approaches  suitable  for  simulating  the  behavior 
of  metals  and  brittle  inorganic  materials,  which  are  the  focus  of  study  here.  For 
additional  information,  the  reader  may  refer  to  extensive  review  articles,  in  particular 
regarding  force  field  models  [52,  53], 

2.1.2  Force  fields 

All-atom  force  fields  are  used  in  molecular  dynamics  simulations  of  metals  and  inor¬ 
ganic  materials  at  the  nanoscale.  Some  classes  of  force-fields  can  be  tuned  and  made 
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applicable  to  a  wide  range  of  materials  whereas  others  exist  for  only  specific  elements 
and  compounds.  For  the  sake  of  brevity,  we  will  only  cover  the  atomistic  force-fields 
used  in  this  work,  i.e.,  the  two-body  Lennard- Jones  (L,J)  and  Morse  potentials,  and, 
the  multi-body  Embedded  Atom  Method  (EAM)  and  Reactive  Force  Field  (RcaxFF) 
potentials. 


2. 1.2.1  Two-body  potentials — LJ,  Morse  and  harmonic 

The  simplest  atom-atom  interactions  are  ones  in  which  the  potential  energy  only 
depends  on  the  distance  between  two  particles.  The  total  energy  is  then  given  by, 

^  N  N 

Utotai = 2  yi  yi  r'<;  ■  •  (2.9) 

*A=i  j=1 

where  p(ry)  is  the  distance  between  particles  i  and  j.  Two-bodv  (pair)  potentials 
must  capture  attraction  at  far  distance  leading  to  creation  of  a  bond,  and  repulsion 
at  very  close  distances  arising  from  the  quantum-mechanical  Pauli  exclusion  principle. 

The  Lennard-Jones  pair  potential  can  be  used  to  model  realistically  the  behavior 
of  noble  gases  (e.g.  argon,  neon).  The  form  of  the  energy  functional  is: 


0(nj)  =  ho 


(2.10) 


where  c0  is  a  measure  of  the  energy  minimum  of  the  pair  potential,  and  a  is  a  measure 
of  the  equilibrium  distance  between  two  atoms,  where  the  force  is  zero.  A  crystal  made 
up  of  LJ  interactions  allows  for  the  formation  and  existence  of  defects,  dislocations 
and  vacancies.  Though  very  rudimentary,  the  LJ  potential  can  also  be  used  as  the 
simplest  model  for  metals  in  some  situations,  and  can  be  fitted  to  the  elastic  constants 
and  lattice  spacing  of  a  metal  crystal.  However,  the  model  has  shortcomings  with 
respect  to  stacking  fault  energies  and  anisotropic  elasticity  of  metals.  Other  simple 
models  for  metals  are  the  Morse  potential  and  the  harmonic  potential. 

The  Morse  energy  functional  is  defined  as, 


<t>(rvj)  =  D  [1  -  exP  (~P  (A,  -  A))]2  ■  (2.11) 
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where  r0  stands  for  the  nearest  neighbor  lattice  spacing  and  D  and  (3  are  additional 
fitting  parameters.  A  fit  of  this  potential  for  different  metals,  can  be  found,  for 
instance  in  [54],  The  Morse  potential  allows  greater  freedom  in  fitting  to  experimental 
properties  than  the  L.l  because  of  the  higher  number  of  parameters. 

In  some  cases,  it  is  advantageous  to  linearize  potentials  around  their  equilibrium 
position;  this  leads  to  the  harmonic  bond  potential: 

1  2 

(j){rl:])  =  a0  +  -k  (rVJ  -  roy  .  (2.12) 

where  k  is  a  spring  constant.  r0  is  the  equilibrium  spacing  between  atoms,  and  a0  is  a 
constant  parameter.  The  harmonic  bond  potential  can  be  augmented  with  harmonic 
3-body  angle  and  4-body  dihedral  potentials.  The  major  strength  of  harmonic  po¬ 
tentials  is  their  computational  simplicity,  and  thus,  the  ability  to  model  systems  of 
large  sizes  of  micron  length  scale. 

The  L.l,  Morse  and  harmonic  potentials  are  also  used  for  testing  £ model'  materi¬ 
als,  i.t.,  to  obtain  general  qualitative  and  order  of  magnitude  insight  into  behavior  of 
different  classes  of  materials.  The  low  number  of  parameters  and  their  relative  mag¬ 
nitude'  can  be  varied  easily  to  study  the  impact  of  atomistic  parameters  on  large-scale 
materials  behavior. 

2. 1.2. 2  EAM  potential 

One  of  the  most  widely  used  force-fields  for  metals  is  the  Embedded  Atom  Method 
(EAM)  and  its  modifications  [55.  53.  56,  57].  The  EAM  is  a  multi-body  force-field 
that  takes  into  account  not  only  interactions  between  pairs  of  atoms,  but;  also  atoms 
with  their  entire  surrounding  neighborhood. 

An  EAM  potential  for  metals  is  typically  given  in  the  form, 

£tot  =  Fi  ( Ph.i )  +  \  foj  (R-ij)  ■  ^2  ^3 ^ 

Ph.i  —  Ej(yp  p]  (Rij)  ■ 

where  Etot  is  the  total  energy  of  the  system.  phti  is  the  density  contribution  at 


atom  i  due  to  remaining  atoms  of  the  system,  F\  (p)  is  the  energy  to  embed  atom 
i  in  the  density  p,  (pVJ  (R)  is  the  pair-pair  interaction  between  atoms  separated  by  a 
distance  R.  The  electron  density  depends  on  the  local  environment  of  the  atom  i,  and 
is  captured  here  by  the  distribution  of  other  metal  atoms  around  the  atom  i  within  a 
cutoff  distance.  It  is  usually  captured  by  a  two-bodv  density  measure  that  describes 
how  the  electron  density  changes  as  a  function  of  the  distance  between  two  atoms. 
The  embedding  function  F  describes  how  the  energy  of  the  atom  depends  on  the  local 
electron  density.  Several  different  analytical  forms  are  used  to  represent  these  two 
functions  for  different  metals.  The  other  half  of  the  potential  is  the  two-body  term  0 
that  captures  the  basic  attraction  or  repulsion  of  two  atoms. 

EAM  potentials  allow  a  much  better  representation  of  the  anisotropy  of  elastic 
properties  and  dislocation  properties  than  two-body  pair  potentials.  They  have  been 
used  successfully  in  modeling  several  FCC  metals  such  as  silver,  gold,  copper,  nickel, 
platinum  and  their  alloys  [53,  58].  They  have  been  widely  used  to  unlock  mechanisms 
and  predict  properties  in  nanoscale  metals  size  regime,  e.g.  in  predicting  strengths 
and  phase  transformations  in  metal  nanowires,  and  crystal  size-strength  relations  in 
nanocrystalline  metals.  They  are,  however,  incapable  of  modeling  effects  of  directional 
bonding,  which  is  important  in  metals  with  some  covalent  character.  To  address 
directional  bonding  in  metals,  modified  EAM  methods  (MEAM)  have  been  proposed 
for  materials  such  as  aluminum,  iron  and  several  others  [59,  60]. 

2. 1.2.3  Reactive  force  field  —  ReaxFF 

Reactive  force  fields  [61,  62,  63]  represent  a  strategy  to  overcome  some  of  the  lim¬ 
its  of  classical  force  fields,  in  particular  their  inability  to  model  chemical  reactions. 
Capturing  the  response  of  chemical  bonds  far  from  equilibrium  turns  out  to  be  very 
important  for  modeling  mechanical  response  in  some  materials. 

Reactive  potentials  start  off  with  a  bond-length  to  bond-order  relationship.  This 
is  used  to  obtain  smooth  transitions  between  different  bond  types,  including  single, 
double  and  triple  bonds.  All  connectivity-dependent  interactions  such  as  valence 
and  torsion  angles  are  formulated  to  be  bond-order  dependent.  Energy  contributions 
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disappear  smoothly  upon  bond  dissociation,  so  that  no  energy  or  force  discontinuities 
appear  during  reactions.  They  also  feature  shielded  nonbonded  interactions  such  as 
Van  der  Waals  and  Coulomb  interactions,  without  discrete  cutoff  distances  to  ensure 
smoothness.  The  method  uses  a  geometry-dependent  charge  equilibration  scheme 
(QEq)  [64]  that  accounts  for  polarization  effects  and  redistribution  of  partial  atomic 
charges  as  a  cluster  of  atoms  changes  its  shape. 


metastable 

configurations 

(BO)f 
BO  =  f(r) 


reactive- nonreactive 
small-strain  approximation 


Figure  2-2:  Concept  of  developing  reactive  potentials  [42].  In  a  nonreactive  potential, 
the  potential  expression  does  not  change  during  deformation  of  the  molecule.  The  reactive 
potential  is  based  on  the  idea  that  the  potential  parameters,  e.g.  the  tangent  spring  con¬ 
stant  k,  depends  on  the  local  atomic  environment.  In  this  example,  the  bond  order  (BO) 
modulates  both  the  tangent  spring  constant  and  the  equilibrium  distance  rg.  This  approach 
enables  one  to  express  a  continuous  energy  landscape  with  several  mctastable  states.  The 
bond  order  can  be  directly  related  to  the  bond  distance  between  atoms. 


A  characteristic  feature  of  the  ReaxFF  force  fields  is  that  all  parameters  of  the 
force  field  are  derived  from  fitting  against  quantum  mechanical  (QM)  data  (density 
functional  theory  (DFT)  calculations).  This  process  is  referred  to  as  force  field  train¬ 
ing  and  the  group  of  QM  experiments  is  called  the  training  set.  The  basic  concept  is 
to  obtain  a  close-as-possible  fit  to  properties  of  several  reactions  and  deformations  of 
groups  of  atoms  from  DFT  calculations  and  ReaxFF  results.  Typically  these  prop¬ 
erties  include  elastic  properties  and  equations  of  state,  surface  energies,  dissociation 
energies,  and  energy  landscapes.  The  absence  of  any  empirical  fitting  to  experimental 
data  makes  this  a  completely  QM-based  force  held. 
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Several  flavors  of  the  ReaxFF  potential  have  been  developed.  The  reactive  poten¬ 
tials.  originally  only  developed  for  hydrocarbons,  have  been  extended  recently  to  cover 
a  wide  range  of  materials,  including  metals,  semiconductors  and  organic  chemistry  in 
biological  systems  such  ms  proteins  [65,  66,  61,  63,  67].  Among  other  applications,  they 
have  been  used  successfully  to  reproduce  properties  ol  several  carbon  nano-systems 
such  as  nanotubes,  buckyballs  and  graphene.  They  have  also  been  used  successfully 
in  reproducing  silicon  oxidation  and  fracture  in  silicon  and  silica  [68,  69.  70,  71]. 

However,  due  to  the  increased  complexity*  of  the  force  field  expressions  and  the 
cost  of  the  charge  equilibration  steps,  the  ReaxFF  simulations  are  usually  between  50 
to  100  times  more  expensive  than  non-rcactive  force  fields,  in  terms  of  computational 
power.  Still  they  are  several  orders  of  magnitude  faster,  and  can  model  systems  much 
larger  than  DFT  calculations,  which  can  also  be  used  to  model  bond  rupture. 

2.1.3  Scaling  and  computational  issues 

All-atom  simulations  are  an  expensive  proposition  as  far  as  computational  resources 
and  time  are  concerned.  Large  number  of  atoms,  and  complexities  of  force-held  ex¬ 
pressions  often  lead  to  the  requirement  for  parallel  computing  techniques  and  the  use 
of  supercomputers  with  several  computing  nodes.  The  state-of-the-art  in  terms  of 
system  size  modeled  have  been  several-billion  atom  simulations  [72,  73],  e.g.  shock 
loading  response  of  crystals  [74]  and  reactions  in  energetic  materials  [75],  Most  of 
these  calculations  have  been  carried  out,  on  parallelized  codes  on  supercomputers 
with  >  1,000  computing  nodes.  The  increase  in  computational  speed,  (now  en  1015 
floating  point  operations  per  second  (PFLOPS)),  over  the  last  few  decades  has  defi¬ 
nitely  helped  in  the  success  of  molecular  dynamics  as  a  widely  used  tool.  However, 
this  size  is  still  a  minuscule  part  of  a  macroscopic  sample  with  approximately  1023 
atoms  (corresponding  to  1  mole).  This  size  is  clearly  not  achievable  today,  and  even  if 
possible,  would  provide  huge  challenges  for  data  handling,  storage,  analysis  and  visu¬ 
alization.  However,  for  many  material  properties  it  is  not  necessary  to  consider  1023 
atoms.  Many  thermodynamic  and  mechanical  properties  can  be  captured  in  systems 
of  thousands  of  atoms  or  much  less. 
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Another  scaling  issue  of  molecular  dynamics  is  of  time-scaling.  The  unit  computa¬ 
tional  step  in  time  is  of  the  order  of  a  femtosecond  in  most  MD  simulations,  to  ensure 
system  energy  conservation  and  the  capture  individual  atomic  vibration  events.  This 
implies  a  MD  simulation  run  over  a  few  days  will  capture  a  system  evolution  over  a 
few  nanoseconds  only.  However,  in  many  applications,  such  femtosecond  accuracy  is 
not  helpful,  as  significant  changes  in  the  system  are  rare  events,  and  a  lot  of  com¬ 
putational  time  is  wasted  following  atomic  vibrations.  This  also  requires  the  use  of 
artifices  to  speed  up  rare  event  encounters,  by  imposing  unrealistically  high  loading 
rates  on  the  system. 

Some  of  these  problems  can  be  overcome  by  methods  for  modeling  across  length 
and  time  scales,  called  multiscale  modeling.  The  core  philosophy  here  is  to  use  combi¬ 
nations  of  quantum,  atomistic  or  continuum  methods  to  tackle  problems  of  large  scale, 
while  keeping  computational  costs  tractable.  An  excellent  example  is  the  problem  of 
chemical  complexity  (e.g.  reactions  at  a  crack  tip)  in  the  deformation  of  a,  material. 
The  reaction  may  be  highly  localized  to  the  crack  tip  region,  requiring  quantum  or 
reactive  force  held  usage  in  that  region,  but  the  effect  on  stress  fields  will  be  large 
length  scale,  requiring  the  entire  system  to  consist  of  thousands  of  atoms  or  more. 


2.2  Multiscale  modeling 

The  properties  of  materials  are  often  determined  by  highly  localized  phenomena  that 
influence  the  material  at  wider  spatial  and  temporal  scales.  Materials  whose  proper¬ 
ties  are  affected  by  such  localized  processes  occurring  at  specific  spatial  and  temporal 
scales  can  be  handled  efficiently  with  multiscale  methods,  based  on  the  idea  of  de¬ 
composing  the  computational  domain  to  reflect  requirements  for  spatial  variations  of 
computational  accuracy,  depending  on  specific  mechanisms  that  occur. 

The  deformation  and  fracture  mechanics  of  materials  is  one  such  area  where  pro¬ 
cesses  occur  over  a  wide  range  of  length  scales,  while  processes  can  be  highly  lo¬ 
calized  at  defects  such  as  grain  boundaries  or  crack  tips.  Chemistry  and  mechanics 
are  typically  considered  independently  for  material  deformation  studies  in  large-scale 
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simulations  using  noil-reactive  force  fields.  However,  the  details  of  bond  breaking  and 
formation  have  been  shown  to  have  significant  influence  on  the  macroscopic  fracture 
mechanics  of  materials,  and  can  not  be  neglected  in  order  to  obtain  predictive  models 
of  the  material  behavior  under  extreme  conditions  and  harsh  environments. 

Multiscale  methods  for  deformation  and  stress  analysis  in  materials  can  be  broadly 
categorized  into  two  types,  depending  on  information  passing  between  different  length 
scales:  hierarchical  and  concurrent  [76].  Hierarchical  schemes  include  MD  meth¬ 
ods  whose  underlying  potentials  are  derived  from  ab-initio  quantum  calculations  e.g. 
ReaxFF  [61,  62],  and  many  EAM  [60]  potentials.  Some  of  these  methods  use  phe¬ 
nomenological  parameters  derived  from  smaller  length  scale  behavior,  used  to  charac¬ 
terize  the  coarser  length  scale.  The  primary  challenge  for  these  methods  is  the  need 
for  complete  knowledge  of  all  relevant  mechanisms  and  processes  at  the  lower  scale. 
Moreover,  the  lowest  length  scale  method  in  these  schemes  has  to  have  no  or  minimal 
empirical  character  like  Density  Functional  Theory  (DFT)  [77.  78.  79]  or  quantum 
chemistry  [80]  methods. 

Concurrent  schemes  link  methods  across  length  scales  simultaneously.  They  com¬ 
municate  through  handshaking  approaches  -  often  featuring  regions  where  two  distinct 
methods  are  used  to  implement  the  transition  from  one  method  to  another.  There 
are  many  systems  and  processes  where  concurrent  modeling  across  scales  is  useful  in 
analysis,  for  example,  in  interface-controlled  materials  phenomena  like  grain  bound¬ 
ary  diffusion,  crack  propagation,  material  embrittlement  or  thin  film  adhesion.  A 
major  challenge  is  to  obtain  smooth  coupling  of  the  disparate  simulation  methods, 
with  no  ghost  forces  or  barriers  at  the  method  interfaces,  to  provide  a  smooth  en¬ 
ergy  landscape.  Several  concurrent  multiscale  schemes  have  been  developed,  among 
them  the  MAAD  method  (microscopic,  atomistic  and  ab-initio  dynamics)  [81.  82] 
for  studying  fracture  in  silicon,  the  quasi-continuum  method  [83,  84]  for  coupling 
atomistic  and  finite-elements,  and  the  CADD  method  (coupled  atomistic  and  dis¬ 
crete  dislocation  method)  [85].  A  recent  scheme  developed  is  a  hybrid  ReaxFF- EAM 
coupling,  developed  within  the  Computational  Materials  Design  Framework  [86,  87], 
a  set  of  computational  tools  that  allows  for  easy  integration  of  simulation  methods 
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across  various  length  scales. 

2.3  Link  to  continuum  state  variables 

It,  is  often  useful  to  be  able  to  derive  continuum  thermodynamic  and  mechanical  state 
variables  from  atomistic  simulations.  Temperature  and  pressure  of  a  system  can  be 
easily  measured  as  statistical  quantities  over  a  certain  time  length  of  a  simulation  [52] . 
Below  are  outlined  how  stress  and  strain  are  measured  from  atomistic  simulations. 

2.3.1  Stress 

The  challenge  in  defining  an  atomistic  stress  tensor  is  relating  it  to  a  continuum  stress 
valid  at  every  point  in  space  is  the  discreteness  of  an  atomistic  materials  simulation. 
The  virial  stress  is  defined  as. 

aO  \7  \  ma%.iVa.j  +  -  (hi  j  -  (2-14) 

V  a  a.f3,a^i3  } 

where  vaj  is  the  velocity  of  atom  a  in  the  direction  i,  rna  is  its  mass,  ra8<i  is  the 
distance  vector  from  atom  a  to  atom  0  in  the  direction  i  and  fa8  •  is  the  force  exerted 
on  atom  a  by  atom  0  in  the  direction  j. 

The  first  term  on  the  right  is  the  kinetic  contribution,  and  the  second  the  contribu¬ 
tion  from  forces  between  the  atoms.  The  virial  stress  needs  to  be  averaged  over  space 
and  time  to  converge  to  the  Cauchy  stress  tensor  of  continuum  mechanics  [88,  89]. 


2.3.2  Strain 

Local  atomistic  atomic  strains  can  be  defined  in  terms  of  the  affine  transformations 
that  transform  the  neighborhood  of  an  atom  from  the  unstrained  configuration  to  the 
strained  state. 

To  fit  a  transformation  in  a  least  squares  sense,  we  seek  the  transformation  matrix 
Jq  at  atom  a  that  minimizes  J2seNa  |  JcTq/?  —  ra/9|"  where  0  runs  over  all  the  neighbors 
of  atom  a  and  rCj  is  the  undeformed  distance  vector  from  atoms  a  to  0  and  ra3  is 
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the  deformed  distance  vector  (both  written  in  column  form  here).  The  solution  of 
this  least  squares  minimization  is. 

U  =  Vi„W 'lpha,  wliereV0  =  £  r»3r»|,  W,  =  £  rjf  J,  (2.15) 

PeNa  fteN, 

where  the  superscript  T  denotes  the  transpose  of  the  matrix.  The  Lagrangian  strain 
matrix  can  then  be  calculated  as 

Va  =  ^(Jjja-I):  (2.16) 

where  I  is  the  identity  matrix  [90]. 


2.4  Visualization  and  data  analysis 

Molecular  dynamics  simulations  produce  large  quantities  of  data,  sometimes  running 
into  several  terabytes.  This  is  because  they  usually  save  snapshots  of  the  system 
state  as  the  simulation  proceeds,  storing  atom  positions,  velocities,  energies,  forces, 
stresses,  and  system  temperatures  and  pressures.  It.  is  necessary  to  have  scripts  and 
programs  that  can  post-process  such  huge  data  sets  to  filter  out  useful  information 
and  display  it  in  a  user-friendly  fashion.  For  example,  a,  billion-atom  simulation  of  a 
single  crystal  under  shock  loading  can  provide  a  very  useful  visualization  if  only  the 
resulting  defects  in  the  crystal  arc  visualized  and  not  the  entire  system. 

For  the  analysis  of  atomistic  simulations  for  mechanical  properties,  measures  like 
strain,  stress  or  potential  energy  of  atoms  arc  important  quantities  that  provide  in¬ 
sight  with  respect  to  continuum  mechanics  theories.  However,  it  is  often  also  useful  to 
post-process  the  data  and  derive  new  quantities  providing  information  about  defect 
structures  within  condensed-matter  systems. 

Next  we  discuss  a  few  examples  for  the  analysis  of  crystal  defects  in  metals  and 
inorganic  crystals  that  will  be  useful  in  the  later  chapters. 
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2.4.1  Energy  method 

An  easy  way  to  ‘see’  in  the  interior  of  solids  and  find  defects  is  to  use  the  energy 
method.  Here  only  atoms  with  potential  energy  greater  than  a  critical  energy  ocr 
above  the  bulk  energy  <pb  are  visualized.  This  method  is  very  effective  for  plotting 
dislocations,  microcracks,  surfaces,  grain  boundaries,  voids  and  other  sites  of  high 
energy.  The  method,  however  faces  difficulties  in  visualization  for  defects  with  very 
low  energies  of  formation,  or  systems  at  high  temperatures,  when  the  thermal  energies 
of  ordered  bulk  atoms  become  comparable  to  certain  defect  energies.  It  is  very  useful 
in  crack  propagation  simulations,  as  it  can  quickly  identify  the  position  of  a  crack  tip 
and  formation  of  new  microcracks/  defects  in  the  vicinity  of  an  existing  crack. 

2.4.2  Slip  vector  analysis 

The  slip  vector  analysis  [91]  provides  visualization  of  slip  planes,  dislocation  cores 
and  stacking  faults  in  crystalline  materials.  It  provides  the  exact  Burgers  vector 
and  direction  of  slip  and  can  also  visualize  incipient  slip  (before  a  dislocation  core  is 
completely  formed).  The  slip  vector  of  an  atom  a  is  defined  as 

's  ay/? 

where  ns  is  the  number  of  slipped  atoms,  sf'3  is  the  vector  difference  of  positions  of 
atoms  a  and  ,8  at  the  current  configuration,  and  is  the  vector  difference  of  these 
atomic  positions  at  no  mechanical  deformation. 

The  slip  vector  analysis  was  originally  used  in  visualizing  defect  structures  formed 
during  the  nanoindentation  of  metals.  It  has  also  proved  very  useful  in  simulations 
of  ductile  fracture. 

2.4.3  Centrosymmetry  parameter 

Centrosymmetry  parameter  [92]  measurement  is  used  for  crystalline  structures  with 
a  center  of  symmetry.  It  uses  the  fact  that  centrosymmetric  crystals  remain  cen- 
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trosymmetric  after  homogeneous  deformation.  Each  atom  has  pairs  of  equally  oppo¬ 
sitely  placed  neighbors  in  a  centrosymmetric  crystal.  This  rule  breaks  down  when 
a  defect  is  in  the  neighborhood  of  the  atom  under  consideration.  The  method  is 
particularly  useful  to  distinguish  different  types  of  defects,  and  to  display  stacking 
faults  (which  are  hard  to  observe  using  the  energy  method).  For  example,  defining 
the  centrosymmctry  parameter  in  a  face-centered  cubic  (FCC)  crystal  as. 


=  £ 


m 


^2  Tk-j  +  hy+e 


k—  i 


(2.18) 


where  rkj  is  the  kth  component  of  the  bond  vector  of  atom  i  with  its  neighbor  atom 
j.  and  rkj+ 6  is  the  same  quantity  with  respect  to  the  opposite  neighbor  in  the  FCC 
crystal.  The  method  can  display  defects  at  high  temperatures  also,  where  the  energy 
analysis  starts  to  fail  due  to  the  thermal  fluctuations  of  the  atoms.  It  has  been  used 
successfully  in  many  simulations  of  FCC  metals  to  display  stacking  faults  and  partial 
and  full  dislocation  cores. 


2.4.4  Common  neighbor  analysis 

The  common  neighbor  analysis  [93]  is  also  used  for  crystalline  structures,  to  classify 
atoms  according  to  their  local  crystalline  structure  (FCC,  BCC.  HCP  etc.).  This 
analysis  method  is  based  on  a  nearest-neighbor  graph,  i.e..  a  mathematical  graph 
structure  consisting  of  edges  that  connect  nearest  neighbor  atoms.  Which  atoms  are 
considered  nearest  neighbors  is  determined  by  a  user-defined  cutoff  distance.  Each 
crystal  structure  exhibits  a  characteristic  local  topology  of  the  nearest-neighbor  graph, 
and  each  atom  is  the  classified  as  belonging  to  a  certain  local  crystal  structure.  It  is 
again  very  useful  for  visualizing  stacking  faults  in  FCC  crystals  as  the  atoms  lying 
within  the  fault  have  a  local  HCP  structure. 

A  wide  variety  of  scripts  in  Python,  C  and  FORTRAN  have  been  written  in  this 
thesis  lor  analyzing  data  from  MD  simulations.  A  wide  range  of  functions  available 
in  MATLAB  have  also  been  used  in  this  thesis  to  use  raw  data  from  simulations  to 
extract  material  properties,  displacement  profdes  and  strain  distributions.  All  scripts 
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used  in  this  thesis  are  provided  in  Appendix  A  with  appropriate  descriptions  within 
the  script. 

2.4.5  Visualization  programs 

Visualization  plays  a  crucial  role  in  the  analysis  of  MD  simulation  results,  as  the  raw 
data  represents  a  collection  of  positions,  velocities  and  accelerations  as  a  function  of 
time  [42],  In  particular,  structural  defects  and  patterns  of  cooperative  atom  motion 
are  difficult  to  analyze.  To  address  this  point,  many  visualization  tools  have  been 
developed  for  displaying  atoms,  molecules  and  assemblies  such  as  nanostructures  and 
bulk  systems.  A  rather  versatile,  powerful  and  widely  used  visualization  tool  is  the 
Visual  Molecular  Dynamics  (VMD)  program  [94],  VMD  enables  one  to  render  com¬ 
plex  molecular  geometries  using  particular  coloring  schemes.  VMD  however  suffers 
from  a  drawback  of  being  unable  to  process  large  systems  with  >  100.  000  atoms. 
Atomeye  [95]  is  another  visualization  program  used  extensively  in  this  thesis  for  vi¬ 
sualizing  large  systems.  Aimed  at  crystalline  materials,  it  possesses  several  tools  for 
the  analyses  of  crystals  and  crystalline  disorder.  Very  large  systems  possessing  0(10 
million  atoms/particles)  are  visualized  using  OpenDX  [96]  and  Ovito  [97].  These 
visualizations  are  often  the  key  to  understanding  complex  correlated  dynamical  pro¬ 
cesses  and  mechanisms  in  analyzing  the  motion  of  atoms,  and  they  represent  a  filter 
to  make  useful  information  visible  and  accessible  for  interpretation. 
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Chapter  3 


Strength  enhancement  through 
bone-inspired  metal-matrix 
nanocomposite  design 


Bone,  a  natural  composite  material,  comprises  of  particularly  ‘weak’  constituent 
phases,  soft  protein  matrix  and  hard  hydroxyapatite  mineral  platelets.  These  con¬ 
stituents  possess  very  poor  mechanical  properties  to  be  used  as  stand-alone  structural 
materials.  Hydroxyapatite  is  as  brittle  as  commercial  ceramics,  and  fractures  catas¬ 
trophically,  whereas  proteins  have  the  same  order  of  stiffness  of  soft  polymers  such 
as  polythene.  In  fact,  the  nanostructural  makeup  of  bone  has  been  shown  to  be  of 
crucial  importance  for  its  superior  mechanical  properties  over  its  constituent  phases, 
providing  high  strength  at  high  stiffness  (see  Chapter  1,  section  1.3).  This  raises, 
as  discussed  in  Chapter  1,  interesting  questions  about  the  viability  of  the  transfer 
of  similar  mechanical  property  enhancement  strategies  to  engineering  materials  and 
structures.  In  particular,  the  application  to  the  design  of  metal-matrix  composite 
materials  seems  promising  due  to  the  availability  of  corresponding  ‘soft’  and  ‘hard’ 
metal  constituents.  The  scope  and  viability  of  the  application  of  bone  nanostructure 
to  functional  metal-matrix  composites,  however,  remains  unresolved.  In  this  chapter, 
a  novel  class  of  biomimetic  nanocomposites  inspired  by  the  structural  motif  found  in 
bone  and  nacre  is  proposed,  and  used  to  formulate  structural  composites  for  meclian- 
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ical  loading  applications.  We  then  report  a  series  of  computational  experiments  using 
molecular  dynamics  simulation  to  study  the  yield  response  of  the  metal-matrix  com¬ 
posite  subject  to  tensile  loading,  and  study  the  changes  in  mechanisms  and  properties 
that  arise. 

3.1  Structure  inspiration  from  nanostructure  of  bone 

Metal-based  nanocomposites  provide  a  great  potential  for  applications  in  high  hard¬ 
ness  and  toughness  material  design.  Potential  applications  include  coatings  for  fric¬ 
tion  and  wear-resistant  cutting  tools,  shock  impact-dissipating  structures  and  other 
tribological  applications  where  strong  functional  materials  are  the  key  to  initiate  fur¬ 
ther  technological  development  [98,  99,  100,  101].  Recent  advances  in  the  development 
of  nanocomposite  materials  have  suggested  that  a  new  paradigm  of  composite  design 
might  be  to  systematically  engineer  the  nanostructural  arrangement  of  components 
by  designing  their  properties,  interfaces  and  geometry  to  tailor  desired  macroscopic 
functional  properties.  These  efforts  extend  earlier  studies  of  creating  nanomaterials 
out  of  metallic  constituents  (e.g.  nanowires,  thin  films)  towards  bulk  materials  [102]. 
However,  the  optimal  choice  of  nanostructural  arrangement  of  material  constituents 
to  maximize  performance  remains  unknown,  preventing  us  from  systematically  car¬ 
rying  out  a  bottom- up  design  approach. 

A  variety  of  biological  structural  materials  such  as  bone  and  nacre  are  known  to 
feature  a  common  “brick-and-mortar”  structural  motifs  at  the  nanoscale,  composed 
of  material  constituents  with  disparate  properties  [103,  33,  46,  23,  7]  (Figure  3-1). 
These  universal  nanostructures  arc  seen  to  combine  inferior  building  materials,  soft 
protein  and  brittle  calcite  or  hydroxyapatite  crystals  to  obtain  structures  with  high 
strength  and  high  toughness  at  biological  scales  [1,  104.  13].  Their  improved  proper¬ 
ties  have  been  attributed  to  their  hierarchical  structure,  as  well  as  their  fundamental 
structural  organization  of  constituting  elements  at  the  nanoscale  [36]  (see  Chapter 
1).  The  biological  role  of  these  materials  is  strongly  related  to  load  carrying  and 
armor  protection  in  nature.  Based  on  their  intriguing  properties,  these  materials 
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raise  an  important,  question  whether  their  design  strategies  could  provide  directions 
for  conventional  structural  engineering  material  design.  However,  for  materials  de¬ 
velopment.  the  use  of  proteins  is  not  a  viable  option,  because  these  materials  are 
rather  difficult  to  synthesize  and  engineer.  Here  we  propose  an  alternative  approach, 
based  on  using  metal-metal  nanocomposites  that  utilize  the  material  concepts  iden¬ 
tified  from  biological  analogs  as  guiding  principles  in  the  design  process.  However, 
despite  earlier  studies  [30,  105.  106],  the  transferability  of  designs  found  in  biological 
structures  towards  conventional  metal  and  ceramic  based  composites  remains  an  un¬ 
resolved  question,  partly  because  the  fundamental  mechanisms  of  how  structure  and 
properties  are  related  have  not  yet  been  explored.  Specifically,  the  wide  parameter 
space  associated  with  different  platelet  shapes  and  orientations  has  not  been  described 
in  the  literature. 

Mimicking  the  universal  design  motif  found  in  mineralized  protein  tissues  (Figure 
3-1)  using  metals  requires  a  systematic  design  of  a  soft  metal /hard  metal  composite 
structure.  To  address  this  issue,  here  we  report  a  systematic  analysis  of  the  deforma¬ 
tion  mechanisms  and  yield  strength  for  a  broad  variation  of  nanostructural  arrange¬ 
ments,  by  using  full  atomistic  molecular  dynamics  simulations.  This  study  elucidates 
the  different  deformation  mechanisms  under  plasticity  and  effect  of  the  nanostructure 
design  parameters  on  their  yielding  behavior.  We  identify  optimal  design  solutions 
by  providing  specific  length-scales  and  geometries  that  yield  maximum  strength  at 
efficient  material  use. 


3.2  Model  construction 

All  simulation  studied  arc  performed  using  fully  atomistic  simulations,  using  EAM 
potentials  (see  section  2.1.2).  All  nanocomposite  structures  under  study  are  con¬ 
structed  out  of  single  crystals  of  the  matrix  and  platelet  phases.  Face-centered  cubic 
(FCC)  metals  are  chosen  for  both  phases,  due  to  availability  of  well-tested  atomistic 
potentials  [55.  107].  Structures  are  designed  such  that  crystal  orientations  of  both 
matrix  and  platelets  are  A' = [Tl 0]  T=[lll]  Z=[U2\  (for  geometr  y  see  Figure  3-1  (b)). 
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To  replicate  uniaxial  tensile  load,  load  is  applied  in  X  direction,  allowing  system  re¬ 
laxation  in  the  Y  and  Z  directions  to  zero  applied  bulk  normal  stresses  in  these  two 
directions  (Figure  3-  1(b)). 


soft  metal  matrix  ^arc*  meta'  platelets 

x  =  [\\0)  (b) 

Figure  3-1:  Panel  (a)  The  ultrastructure  of  bone  showing  hard  mineral  platelets  2-4  run 
thick  and  up  to  100  nm  long  embedded  in  a  soft  collagen-rich  matrix  (figure  adapted  and 
redrawn  from  reference  [33]).  Panel  (b):  A  geometrically  two-dimensional,  simple  schematic 
model  of  hard/soft  phases  arranged  in  a  nanocomposite,  based  on  the  ultrastructure  of 
bone.  This  structure  is  realized  with  metallic  components  using  a  “soft”  Cu  metal  matrix 
and  modified  “hard”  EAM  model  metal  platelets.  To  create  this  composite  structure,  a 
single  copper  crystal  is  arranged  in  the  [110]  x  [111]  ( X-  Y )  orientation,  and  a  regular  array 
of  rectangular  voids  is  created.  The  platelet  crystals  are  inserted  in  the  voids  in  the  same 
crystallographic  orientation,  ensuring  no  overlapping  atoms  (by  avoiding  any  distance  closer 
than  Cu-Cu  nearest  neighbors).  The  resulting  structures  are  relaxed  to  minimize  the  global 
normal  stresses.  The  plot  also  shows  the  various  parameters  that  determine  the  geometrical 
arrangement  of  such  a  system.  Sec  Eq.  3.2  for  the  relation  between  these  parameters  and 
the  ones  used  in  our  study. 


To  implement  a  fully  atomic-scale  study  of  deformation  in  these  materials,  atom¬ 
istic  interaction  potentials  are  required  for  an  accurate  description  of  deformation  in 
the  soft  matrix  phase,  the  hard  platelet  phase,  and  a  cross-potential  for  interactions 
across  the  matrix-platelet,  interfaces.  The  embedded  atom  method  (EAM)  alloy  po¬ 
tential  [53,  58]  is  chosen  to  model  the  two  phases  and  their  interactions.  The  EAM 
model  has  been  shown  extensively  to  provide  an  accurate  description  of  FCC  met- 


als  and  their  alloys  [56].  The  use  of  the  BAM  alloy  potential  allows  the  freedom  of 
modeling  of  two  distinct  metals  with  quite  different  stiffness  and  strength  properties. 

The  samples  are  subject  to  energy  minimization  to  relax  internal  stresses,  followed 
by  quasi-static  tensile  loading  using  molecular  statics.  System  sizes  chosen  for  simu¬ 
lations  are  approximately  500  A  x  200  A  x  25  A.  with  approximately  200,000  atoms 
in  each  system.  Periodic  boundary  conditions  are  applied  in  all  three  directions,  to 
mimic  large  crystal  systems.  The  system  size  is  large  enough  to  allow  dislocation 
interactions  across  4  layers  of  platelets,  but  not  too  large  to  slow  down  computa¬ 
tion.  Larger  sized  samples  with  periodic  boundary  conditions  show  similar  flow  stress 
values.  Load  application  is  achieved  in  two  steps-  system  global  stress  relaxation, 
followed  by  application  of  uniaxial  tensile  load.  Quasistatic  loading  is  achieved  using 
displacement  boundary  conditions  along  with  an  energy  minimization  scheme  imple¬ 
mented  using  a  micro-convergence  integrator.  Flow  stresses  are  calculated  as  system 
averaged  virial  stresses  [52]  at  large  strains,  where  the  system  stress  fluctuates  around 
constant  values.  The  flow  stress  calculation  is  averaged  over  a  strain  range  from  0.15 
to  0.2,  over  a  range  where  plasticity  has  initiated,  and  the  stress  is  about  constant. 
The  local  stress  distribution  in  the  sample  is  calculated  through  atomic  level  virial 
stresses  [89]. 

3.3  Interatomic  potential  development  and  testing 

The  model  developed  here  is  not  designed  to  represent  a  particular  pair  of  materials. 
Rather,  it  is  developed  to  allow  the  freedom  to  modify  material  properties  in  their 
absolute  values  and  ratios  for  matrix,  platelets  and  interfaces,  in  the  spirit  of  a  model 
material  [108].  The  concept  behind  using  such  an  approach  is  driven  by  the  desire 
to  explore  a  wide  parameter  space  and  to  understand  the  material  behavior  of  a 
wide  range  of  material  properties,  without  focusing  on  a  single,  specific  material. 
This  facilitates  a  computational  engineering  approach  in  which  we  systematically 
investigate  the  sensitivities  of  key  design  parameters  on  the  overall  material  behavior 
to  provide  generic  understanding  of  the  relationship  between  structure  and  properties. 


We  specific  avoid  focusing  on  modeling  a  specific  material.  The  atomistic  model  in  this 
study  is  based  on  earlier  EAM  potentials  developed  for  two  species  alloy  systems  such 
as  Ni-Al.  Cu,  and  Ag  alloys  [53,  109].  The  Baskes-Daw  model  of  an  EAM  potential 
for  face-centered  cubic  (FCC)  metals  is  reproduced  here  from  Chapter  2,  since,  it  will 
be  used  for  the  modified  potential  development  (see  cq.  (3.1)  below).  It  consists  of 
a  pair  interaction  term,  and  an  electron  density  term  which  contributes  through  an 
embedding  energy  term  that  charges  an  energy  penalty  for  deviating  from  an  FCC 
environment.  The  development  of  an  alloy  potential  requires  the  development  of  a 
cross-pair  potential,  and  density  functions  for  atoms  of  type  A  in  an  environment  of 
type  B,  and  vice  versa.  The  total  potential  energy  is  given  by, 

Em  —  Sj  Ei  ( Ph.i )  +  2  Yli  °ij  ( Eij )  .  ^ 

Ph.i  =  Pj  (Rij)  : 

where  Etot  is  the  total  energy  of  the  system,  ph%i  is  the  density  contribution  at 
atom  i  due  to  remaining  atoms  of  the  system.  F2  (p)  is  the  energy  to  embed  atom 
i  in  the  density  p,  Ay  (R)  is  the  pair-pair  interaction  between  atoms  separated  by  a 
distance  R. 

The  matrix  material  is  chosen  as  a  soft  metal,  copper,  modeled  using  the  Mishin 
potential  [107].  A  model  material  of  higher  stiffness  and  higher  strength  is  chosen  as 
the  platelet  phase.  To  achieve  this,  a  modified  copper  EAM  potential  is  used  to  model 
elemental  second  phase  material.  This  provides  the  freedom  to  modify  stiffness  and 
strength  ratios  between  matrix  and  platelets.  The  inter-elemental  potential  terms 
consist  of  the  cross-potential  pair  potential,  which  has  also  been  chosen  as  a  modified 
Mishin  potential.  Use  of  a  modifiable  cross-potential  allows  us  to  vary  the  interfacial 
strength,  and  effect  of  the  same  on  flow  stress.  To  design  the  modified  interatomic 
potentials,  we  selectively  change  the  pair  interaction  of  the  copper  EAM  potential 
while  keeping  the  density  and  embedding  energy  terms  unmodified.  Modification 
of  only  the  pair  term  allows  us  to  modify  the  bulk  modulus,  cohesive  energy,  and 
unstable  and  stable  stacking  fault  energies,  while  maintaining  lattice  parameter  of 
the  FCC  phase  constant.  We  use  the  modified  potentials  for  describing  the  platelets 
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and  interface,  which  yields  the  potential  properties  shown  in  Figure  3-2.  The  ratio  of 
elastic  moduli  of  the  metal  phases,  are  not.  as  dissimilar  as  for  bone  constituents,  due 
to  limitations  of  the  materials  used. 

3.3.1  Interatomic  potential  properties 

The  cross-species  potential  parameters  are  adjusted  so  that  the  matrix-platelet  in¬ 
terface  is  always  weaker  than  both  the  matrix  and  platelets.  In  the  present  study, 
two  different  intcrfacial  strength  levels  are  used,  referred  to  as  “strong5*  and  “weak" 
interface,  defined  in  this  fashion  by  the  adhesion  energy  across  the  (111)  interface. 
The  work  of  adhesion  for  these  two  interfaces  in  the  (111)  orientation  is  compared 
with  several  example  metal/ceramic  interfaces  from  literature  in  Table  3.1  (where 
the  “strong”  interface  resembles  the  surface  energies  of  an  Al/TiC  structure,  and 
the  "weak”  interface  a  Cu/AFCF  structure).  The  freedom  in  variation  of  interfacial 
strength  allows  mechanisms  of  matrix-platelet  decohesion  and  interfacial  sliding  to 
be  activated  and  become  dominant  deformation  mechanisms  at  different  stress  levels. 
The  strength  and  stiffness  of  the  matrix,  platelets,  and  interfaces  are  summarized  in 
Figure  3-2a.  The  stiffness  of  the  platelet  material  is  approximately  two  times  that 
of  the  matrix  material.  This  is  close  to  actual  stiffness  ratios  of  soft  and  stiff  metals 
such  as  Cu/Fe  and  Cu/Pt.  The  generalized  stacking  fault  curves  [110]  have  also  been 
calculated  for  the  matrix,  platelets  and  interfaces,  and  are  shown  in  Figure  3-2b.  The 
large  unstable  stacking  fault  energy  [111]  for  the  platelet  material  shows  the  difficulty 
of  nucleating  dislocations  in  this  material,  thus  making  it  a  “strong”,  less  ductile 
material. 


Interface  type 

W  adhesion  (in  J / HI " ) 

Strong  interface 

2.42 

Weak  interface 

1.12 

Al/TiC 

2.63 

Cu/A1203 

1.09-1.30 

Table  3.1:  Work  of  adhesion  across  a,  “strong”  and  “weak”  interface  and  comparison  to 
actual  metal-ceramic  interfaces. 
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Figure  3-2:  Properties  of  the  interatomic  potentials,  for  the  various  cases  considered 
here.  Panel  (a)  shows  the  theoretical  shear  strength  against  (111)  plane  shearing  in  the 
[112]  direction,  for  matrix  (Mishin  Cu  potential),  platelet,  “strong"  interface  (interfacel) 
and  "weak'’  interface  (interface2)  EAM  potentials,  showing  matrix:platelet  stiffness  ratio 
of  1:2.  Panel  (b)  depicts  generalized  stacking  fault  curves  for  the  same,  showing  high 
dislocation  nucleation  barrier  in  platelet  phase.  Panel  (c)  shows  enthalpies  of  formation  of 
intcrmetallics  and  heats  of  mixing  of  disordered  alloys  between  matrix  (A)  and  platelet  (B) 
materials  at  both  “strong”  and  “weak”  interfaces,  showing  very  high  positive  values  leading 
to  low  solubility  of  one  material  in  another  [106]. 


56 


The  thermodynamic  stability  of  atomically  sharp  interfaces  between  matrix  and 
platelets  as  they  appear  in  these  models  is  also  investigated.  The  heat  of  formation 
has  been  calculated  for  alloys  of  matrix  and  platelets  materials  with  both  types  of 
interfaces  over  the  entire  composition  range.  The  heat  of  formation  of  possible  inter- 
metallic  structures  AB3  and  A3B  (where  A  is  matrix  and  B  is  platelet  material)  arc 
calculated  assuming  a  Ll0  structure,  and  for  AB  assuming  a  Ll2  structure  (Figure 
3-2c)  [57,  112],  The  heats  of  formation  of  the  disordered  alloys  over  the  entire  com¬ 
position  range  are  calculated  from  these  values  using  the  cluster  expansion  method 
(Figure  3-2c)  [113].  The  large  positive  heats  of  formation  for  the  intermet, allies  show 
no  possibility  of  them  forming  at  the  interfaces.  The  positive  and  large  heats  of 
formation  of  the  disordered  alloys  also  hint  at  very  low  solubility  of  one  element  in 
another  and  vice  versa.  Thus  overall  the  implicit  assumption  of  stable  interfaces 
between  matrix  and  embedded  platelets  is  justified. 


3.3.2  Design  Parameters 

The  quasi-two  dimensional  geometrical  arrangement  of  platelets  is  characterized  uniquely 
in  terms  of  five  independent  geometric  parameters  (shown  in  Figure  3- lb).  These  pa¬ 
rameters  are  assuming  rectangle  shaped  platelets  arranged  on  a  regular  two-dimensional 
motif.  Irregular  shaped  platelets  can  also  be  considered  in  terms  of  average  platelet 
sizes  and  distances,  to  be  specified  by  the  five  independent  parameters.  The  relation 
between  different  parameters  is  given  below: 

T  /  _ lx  ly _ 

f  (lx  -\-Wx  )  (ly  ~\~Wy  )  7 

A  -  IJly,  (3.2) 

A  =  l  1 

where  Vf,  lx,  ly.  wx,  wy.  A  and  Ap  are  volume  fraction,  length  and  width  of  platelets, 
axial  and  transverse  spacing  of  platelets,  aspect  ratio,  and  area  of  a  single  platelet. 
The  axial  and  transverse  spacing  between  platelets  are  defined  as  spacing  parallel  and 
perpendicular  to  the  loading  direction  X  respectively.  The  independent  parameters 
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chosen  here  are  Vf,  A,  wx,  wy  and  /,  where  /  is  defined  as  the  stagger  between 
platelets  across  successive  rows  as  a  fraction  of  the  axial  repeat  distance. 

The  choice  of  these  parameters  helps  us  to  answer  questions  as:  are  elongated 
platelets  needed  and  why?  How  far  apart  should  the  platelets  be  spaced?  What  is  the 
effect  of  volume  fraction  of  platelet  phase  on  flow  stress?  Changes  in  other  parameters 
can  be  studied  as  a  combination  of  changes  of  any  of  these  five.  First  order  effects 
of  each  of  these  parameters  are  studied  by  varying  one  particular  parameter  while 
keeping  the  other  four  constant. 

3.4  Atomistic  simulations  under  tensile  deforma¬ 
tion 

3.4.1  Effect  of  geometric  parameters 

Studies  of  effect  on  flow  stress  of  varying  all  geometric  parameters  are  presented  in 
the  subsequent  sections. 

3.4. 1.1  Effect  of  platelet  volume  fraction  Vf 

First  the  platelet  volume  fraction  is  varied  from  «10%  to  «45%,  and  effect  on  flow 
stress  is  measured.  Keeping  the  aspect  ratio  and  spacing  between  platelets  constant, 
the  increase  in  volume  fraction  is  achieved  by  increasing  absolute  size  of  the  platelets. 
There  is  a  slight  decrease  in  flow  stress  as  volume  fraction  is  increased,  but  the  effect 
is  not  very  strong.  The  results  are  shown  in  Figure  3-3. 

3. 4. 1.2  Effect  of  platelet  offset  / 

Platelet  rows  parallel  to  the  axial  direction  have  a  pronounced  offset  from  one  row 
to  the  next  in  structures  such  as  nacre  and  bone  (see  Figure  3- la).  To  explore  the 
significance  of  this  offset,  a  systematic  analysis  using  the  present  model  system  is 
carried  out.  The  offset  in  the  present  experiments  is  measured  as  a  fraction  of  the 
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Figure  3-3:  Effect  of  volume  fraction  on  flow  stress.  Volume  fraction  is  varied  from  ~  10% 
to  45%  and  a  slight  dependence  of  flow  stress  is  seen  within  this  range.  Black  line  depicts 
best  linear  fit  to  data.  Insets  show  nanostructure  at  different  volume  fractions.  The  flow 
stress  of  the  pure  copper  matrix  is  shown  for  comparison  (triangular  marker  on  the  left,  with 
1.800  parts  per  million  platelet,  material,  to  induce  heterogeneous  dislocation  nuclcation). 

repeat  unit  length  in  the  axial  direction,  and  thus,  can  vary  from  0  to  1.  The  flow 
stress  is  measured  as  a  function  of  this  offset  and  results  are  shown  in  Figure  3-4a. 

A  significant  change  of  ~10  times  is  identified  in  the  flow  stress,  with  a  saturation 
observed  at  a  platelet  offset  of  «0.25.  The  load  distribution  within  the  structure  is 
studied  through  the  local  shear  stress  distribution  axy  in  the  regime  of  elastic  de¬ 
formation  (just  prior  to  plasticity),  as  depicted  in  Figure  3-4b.  Significant  effects  of 
the  details  of  shear  stress  distribution  as  a  function  of  nanostructured  geometry  is 
observed,  with  large  homogeneous  domains  of  shear  stress  transfer  found  at  platelet 
offsets  beyond  0.25.  Gao  s  shear- tension  load  sharing  model  [33]  is  seen  to  be  ap¬ 
plicable  to  this  case,  since  at  offsets  of  «0.25  to  «0.75  there  is  sufficient  overlap 
between  platelets  in  subsequent  rows  to  induce  predominant  shear  loading  of  the 
matrix  material  situated  between  platelet  rows.  This  optimal  platelet  overlap  leads 
to  an  effective  distribution  of  load  that  plays  to  the  respective  strengths  of  platelet 
and  matrix  material.  Specifically,  platelets  are  subjected  to  tensile  load  with  very 
small  shear  loading,  and  the  ductile  matrix  material  is  under  shear  loading  with  very 
small  tensile  loading.  These  results  show  a  critical  platelet  offset  beyond  which  the 
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Figure  3-4:  Effect  of  platelet  offset  on  flow  stress,  (a)  shows  large  effect  of  variation  of 
offset  on  flow  stress,  and  a  saturation  of  value  at  w0.25  offset;  insets  show  the  nanostructure 
at  different  offset  amounts  (b)  shows  the  X-  Y  plane  local  shear  stress  distribution  in  the 
nanostructure,  with  dashed-line  bars  indicating  the  location  of  platelet.  The  results  show 
that  effective  load  transfer  between  platelets  through  shear  of  matrix  material  occurs  at 
offsets  larger  than  «0.25  (note  the  symmetry  of  the  problem). 


shear-tension  model  provides  an  accurate  model. 


3. 4. 1.3  Effect  of  aspect  ratio  A 

The  aspect  ratio  of  second  phase  particles  plays  an  important  role  in  load  transfei 
and  anisotropic  properties  in  composites.  The  discontinuous  reinforcement  material 
shape  can  vary  from  a  needle  geometry  (oriented  in  the  direction  of  loading)  to  a  flat 
disk  shape  in  the  other  extreme.  Platelets  in  nacre  and  bone  are  typically  observed 
to  be  of  flat  disk  shape  with  large  aspect  ratios  in  length/thickness  to  width. 

7  T 
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Figure  3-5:  Effect  of  platelet  aspect  ratio  on  flow  stress.  The  results  reveal  that  large 
aspect  ratio  platelets  provide  a  better  mechanical  performance,  as  the  How  stress  saturates 
at  a  maximum  value.  Insets  show  the  nanostructure  at  different  aspect  ratio  values. 


In  the  two-dimensional  geometric  nanostructure,  we  measure  the  aspect  ratio  as 
the  relation  of  length  of  platelets  in  the  axial  direction  to  the  length  in  the  transverse 
direction.  The  aspect  ratio  is  varied  from  ~1  (square  platelets)  to  «20  (highly  elon¬ 
gated  platelets).  The  results  of  the  effect  of  aspect  ratio  on  flow  stress,  while  keeping 
other  four  parameters  fixed,  are  shown  in  Figure  3-5.  For  small  aspect  ratios,  the 
change  of  platelet  ratios  results  in  an  increase  of  the  yield  strength.  A  critical  aspect 
ratio  is  observed  beyond  which  the  flow  stress  is  almost  constant,  occurring  at  an 
aspect,  ratio  of  «4.  An  analysis  of  t  he  tensile  stress  distribution  (results  not  shown) 
reveal  a  significant  change  in  load  redistribution  as  the  aspect  ratio  is  increased,  which 
suggests  that  elongated  platelets  are  to  be  preferred  for  maximal  tensile  load  to  be 
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carried  by  the  platelets. 


3. 4. 1.4  Effect  of  axial  platelet  spacing  wx 

As  the  next  step  we  vary  the  spacing  between  platelets  in  the  same  row.  while  keep¬ 
ing  the  volume  fraction,  platelet  offset,  aspect  ratio  and  transverse  platelet  spacing 
constant.  The  analysis  of  the  flow  stress  as  a  function  of  spacing  shows  a  relatively 
weak  dependency.  However,  despite  the  limited  effect  on  the  flow  stress,  a  change 
in  deformation  mechanism  is  identified.  At  small  spacing,  we  observe  platelet-matrix 
decohesion,  and  partial  dislocation  emission  from  interfacial  cracks  (Figure  3-6).  At 
larger  spacing,  no  decohesion  is  found,  and  dislocations  aic  emitted  bom  legions  ol 
stress  concentration,  such  as  the  corners  of  platelets. 


1  1 
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Figure  3-6:  Effect  of  axial  platelet  spacing  wx  on  the  flow  stress.  Axial  spacing  is  the 
distance  between  platelets  along  the  loading  axis.  Minimal  cftcct  on  flow  sticss  value  is 
seen,  though,  change  in  mechanism  is  observed  as  the  spacing  is  increased.  The  upper 
insets  shows  cent  rosy  mmetry  plots  of  slipped  regions  during  initial  plasticity  with  coloied 
bars  indicating  location  of  platelets,  and  the  lower  insets  show  corresponding  nanostructure. 
At  lower  spacing,  platelet-matrix  decohesion  occurs,  shown  by  the  red  circles,  followed  by 
dislocation  emission  from  the  cracks.  At  higher  spacing,  dislocation  emission  is  seen  from 
areas  of  stress  concentration  such  as  platelet  corners. 
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3.4. 1.5  Effect  of  transverse  platelet  spacing  wy 


The  spacing  transverse  to  the  direction  of  loading  is  varied  in  the  next  step,  while 
keeping  the  other  four  parameters  fixed.  These  analyses  show  a  gradual  decrease  in 
the  flow  stress  as  the  transverse  spacing  is  increased.  A  plot  of  the  parameter 
versus  flow  stress  shows  a  linear  trend,  suggesting  Hall-Petcli  like  behavior  where 
the  flow  stress  increases  by  «50%  [114,  115]  (Figure  3-7).  The  analysis  reveals  that 
the  effective  length-scale  parameter  for  the  Hall-Pctch  like  behavior  is  the  transverse 
spacing.  wy.  This  is  probably  due  to  dislocation  motion  barrier  established  at  the 
matrix-platelet,  interface,  which  effectively  acts  as  a  grain  boundary-like  stiucture. 


1/\Wy  (1/  A) 


Figure  3-7:  Effect  of  transverse  spacing  wy  between  platelets  on  flow  stress.  The  plot 
shows  the  flow  stress  as  a  function  of  1/^,  revealing  a  Hall-Pctch  like  behavior  up  to  the 
strength  limit  (insets  show  nanostructuie  at  dilleient  spacings). 


3. 4. 1.6  Effect  of  intcrfacial  strength 

The  effect  of  matrix-platelet  intcrfacial  strength  on  deformation  mechanisms  and  flow 
stress  is  studied  by  performing  a  self-similar  size  scale  analysis  for  both  the  “weak 
and  “strong”  interfaces  (see  Table  3.1).  The  non-dimensional  parameters,  that  is, 
platelet  offset,  aspect  ratio  and  volume  fraction  are  kept  constant,  while  the  size  of 
the  platelets  Ap  is  scaled  up.  Figure  3-8  shows  how  t  he  flow  stress  varies  with  1/ y/w^ 
for  both  the  interface  types. 
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Figure  3-8:  Effect  of  interfacial  strength  on  flow  stress  by  varying  the  composite  size, 
by  scaling  dimensions  in  same  ratio.  Panel  (a)  displays  effect  of  system  size  on  flow  stress 
for  “strong”  interface,  showing  a  Hall-Petch  like  behavior  up  to  a  strength  limit.  Panel  (b) 
shows  the  effect  of  system  size  on  flow  stress  for  “weak  interface  showing  transition  from 
Hall-Petch  like  to  inverse  Hall-Petch  like  behavior  at  small  scales.  The  circled  region  shows 
the  length  scale  at  which  this  transition  occurs. 
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We  observe  a  significant  qualitative  difference  in  behavior,  depending  on  the  type 
of  interface  between  the  platelet,  and  matrix.  The  structure  with  “weak”  interface 
has  an  optimal  size  at  which  the  flow  stress  is  maximal,  whereas,  the  structure  with 
“strong"  interface  reaches  a  strength  limit  at  lower  sizes.  The  structure  with  “weak” 
interfaces  features  dislocation  sliding,  and  interfacial  debonding  mechanisms  activated 
under  plastic  deformation.  At  small  sizes,  these  become  the  dominant  mechanisms 
due  to  the  increase  in  interfacial  area  in  the  structure.  This  is  reminiscent  of  the 
grain  size  dependence  of  yield  stress  in  nanocrystalline  metals,  and  hence  the  small 
size  regime  is  referred  to  as  inverse  Hall-Petch  behavior  in  Figure  3-8b.  The  “strong” 
interface  structure,  on  the  other  hand,  has  only  bulk  dislocation  mechanisms  for 
plasticity. 

The  study  of  interfacial  strength  suggests  that  the  optimal  size  of  platelets  for 
maximum  flow  stress  strongly  depends  on  the  interfacial  strength.  For  interfaces  of 
high  strength,  which  suppress  interfacial  sliding  altogether,  the  results  suggest  that 
there  exists  a  critical  size  below  which  the  flow  stress  is  size-independent,. 

3.4.2  Sensitivity  analysis  of  design  parameters 

The  sensitivity  of  changes  in  flow  stress  to  changes  in  any  of  the  five  parameters  is 
analyzed  to  determine  the  important  ones  which  greatly  influence  the  yielding  behav¬ 
ior.  The  sensitivities  are  measured  as  rates  of  changes  of  flow  stress  with  parameter 
change,  and  results  are  shown  in  Figure  3-9.  The  figure  illustrates  that  platelet  off¬ 
set  and  transverse  spacing  are  the  most  significant  parameters  that  control  the  flow 
stress.  The  platelet  offset  affects  the  flow  stress  by  changing  the  local  load  distri¬ 
bution,  whereas  transverse  spacing  has  an  affect  through  a  Hall-Petch  grain  size-like 
effect. 

3.5  Strength  saturation  with  size 

We  observe  that  the  interfacial  strength  between  matrix  and  platelets  is  a  key  pa¬ 
rameter  that  determines  the  importance  of  sliding  and  decohesion  as  deformation 
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Figure  3-9:  Sensitivity  of  flow  stress  to  different  design  parameters  as  measured  by 
rate  of  change  of  flow  stress  to  parameter  value  change.  Panel  (a)  shows  sensitivity  limits 
(minimum-to-maximum)  to  non-dimensional  parameters,  that  is,  volume  fraction,  platelet 
offset  and  aspect  ratio.  Panel  (b)  shows  the  sensitivity  limits  of  dimensional  parameters. 
i.e..  axial  spacing  and  transverse  spacing. 
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mechanisms'  at  small  spacing  and  platelet  sizes.  This  leads  to  a  peak  flow  strength  for 
the  nanocomposite  as  a  function  of  platelet  size  (observed  at  a  critical  length  scale  of 
17  A  for  the  "weak”  interface  considered  here).  However,  when  the  interface  is  strong 
enough  to  prevent  sliding  and  platelet  pullout  well  after  the  matrix  has  started  de¬ 
forming  plastically,  we  observe  no  such  peak  flow  strength.  Instead,  we  see  a  strength 
saturation  for  platelet  spacing  below  a  certain  value.  This  raises  interesting  questions 
about  what  happens  to  plasticity  mechanisms  in  the  matrix  under  confined  conditions 
(bound  by  interfaces  with  platelets  on  top  and  bottom).  To  understand  this  strength 
saturation  more  carefully,  we  devise  experiments  to  study  the  effect  of  specimen  size 
confinement  on  dislocation  plasticity  mechanisms  in  the  following  sections. 


3.6  The  size  confinement  effect  on  dislocation  plas¬ 
ticity 

Confinement  effects  have  been  previously  suggested  as  an  explanation  for  the  high 
strength  of  nanostructured  materials.  It  has  been  suggested  that  the  underlying  phys¬ 
ical  reason  is  that  the  nueleation  and  propagation  of  dislocations  becomes  increasingly 
difficult  under  reduction  of  the  material  size.  For  example,  size  effects  in  nanocrys¬ 
talline  materials  [116,  117.  118.  119]  and  nanostructured  ultra-thin  films  [120,  121] 
suggest  fundamental  changes  in  deformation  mechanisms  once  the  dimensions  of  the 
microstructure  approach  nanoscale,  often  characterized  by  the  breakdown  of  complete 
dislocation  mechanisms  [122,  123,  124]  and  dominance  of  interfacial  processes. 

However,  Tip  until  now  a  systematic,  joint  atomistic-nanomechanical  analysis  of 
size  effects  of  ductile  materials  has  not  been  reported.  Earlier  studies  have  been 
carried  out  by  explicitly  including  the  material  nanostructure  (e.g.  polycrystallinc 
metals  [116.  117,  119]).  However,  due  to  the  complexities  of  the  nanostructures  in 
these  materials  these  studies  have  not  yet  enabled  a  systematic  investigation  of  size 
effects  limited  purely  to  dislocation  processes,  since  other  competing  phenomena  (e.g. 
grain  boundary  mechanisms)  have  been  present. 
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Here  we  investigate  the  atomistic  mechanisms  of  plastic  deformation  of  a  ductile 
single  crystal  under  size  confinement  effects  via  a  combination  of  theoretical  and 
molecular  dynamics  analyses.  The  goal  of  this  analysis  is  to  show  that  geometric 
confinement,  of  materials  can  control  the  deformation  mechanism,  once  the  material 
dimensions  reach  a  characteristic  length  scale.  Wo  design  a  simplistic  thin  strip 
single  crystal  model  system  that  enables  us  to  focus  solely  on  size  effect  studies  of 
mechanisms  of  full  and  partial  dislocation  plasticity.  This  is  admittedly  simple  model 
system,  however,  enables  us  to  investigate  fundamental  length  scale  limits  associated 
with  dislocation  mechanisms. 


3.7  Theoretical  analysis  of  the  size  confinement  ef¬ 
fect 


Our  theoretical  analysis  is  based  on  the  Rice-Peierls  model-]  11 1]  that  predicts  a  critical 
energy  release  rate  for  nuclcation  of  leading  and  trailing  partial  dislocations  from  a 
crack  tip  as  a  function  of  the  unstable  stacking  fault  energy  7us  and  the  stacking  fault 
energy  -ty  k-  We  consider  a  thin  strip  geometry  with  width  £  and  a  semi-infinite  crack 
of  length  a  »  £  on  a  single  {111}  slip  plane,  as  shown  in  Figure  3-10.  The  choice  of 
this  model  is  motivated  by  the  fact  that  this  geometry  is  accessible  to  both  theoretical 
and  molecular  dynamics  studies.  In  this  geometry,  the  semi-infinite  crack  plane  and 
slip  plane  in  an  elastically  isotropic  material  coincide,  and  a  dissociated  dislocation 
will  move  out  from  the  crack  tip  in  the  slip  plane  under  pure  mode  II  loading.  For  this 
geometry,  the  critical  energy  release  rates  for  nucleation  of  the  leading  and  trailing 
partial  dislocations  from  the  crack  tip  are  [111]: 

^jli  =  (l  +  (1  —  v)  tan'  pus,  (3.3) 
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where  v  is  Poisson's  ratio,  p  is  the  shear  modulus,  bA  is  the  partial  dislocation 
Burgers  vector  length,  and  rA  denotes  the  separation  distance  between  the  two  partial 
dislocations,  and  oA,  <Pb  are  the  angles  between  the  part.ials  Burgers  vectors  and  the 
X-axis  in  the  slip  plane  [111]. 


I 


crack  tip 


Figure  3-10:  Configuration  and  atomistic  interactions  of  the  thin  strip  geometry.  Crack 
length  from  left  end  of  sample  to  crack  tip  is  300  A.  Morse  interactions  are  defined  across 
the  slip  plane/  crack  plane,  i.c..  for  pairs  of  atoms  in  regions  1-4.  2-3  and  2-4;  there  arc  no 
interactions  across  the  slip  plane  for  regions  1-3  (signifying  the  crack  region);  and  all  other 
interatomic  interactions  are  harmonic  pair  potentials. 
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We  consider  this  model  crack-slip  plane  system  bounded  within  a  slab  of  height  £ 
(see  Figure  3-10).  The  energy  release  rate  G  for  this  thin  slab  geometry  is  [125]: 


(3.6) 


Setting  G  equal  to  Gc[Iltl  (that  is.  we  require  that  GfIiti  =  G)  leads  to  a  critical 
load  T„crit  = 


\i  ~  l,  2). 

Since  rGrit  ~  if y/£,  the  strength  of  the  crystal  increases  as  its  size  is  reduced. 
However,  the  strength  of  the  single  crystal  must  be  limited  by  the  theoretical  shear 
strength  rlh,  the  strength  limit  associated  with  homogeneously  shearing  the  crystal. 
Thus  at  small  dimensions,  the  crystal  must  fail  by  homogeneous  shear  across  the  crack 
plane,  representing  a  breakdown  of  the  ability  to  fail  under  dislocation  nuclcation. 

We  now  consider  an  FCC  crystal  with  a  [110]  (/'-direction  and  a  (111)  crack  plane. 
This  implies  that  the  two  partials  to  emerge  in  the  (111)  slip  plane  will  be  at  ±30 
degrees  orientation  relative  to  the  x-axis  in  the  x-z  plane.  By  setting  r  =  rth.  we 
obtain  two  critical  length  scales, 


i?"  =  ^GJt/Tpcc  and  (3.7) 

We  note  that  after  the  first  partial  or  uniform  shear  event,  the  slip  plane  contains 
an  HCP  stacking  fault  region,  with  t^cp  ^  TtPCC  (Figure  3-11  (b)). 

This  model  provides  important  predictions:  At  large  crystal  dimensions  £  >  ■ 

shear  deformation  is  mediated  by  nucleation  of  full  dislocations,  representing  the 
conventional  regime  of  plasticity.  However,  below  £  <  £Wb  deformation  proceeds 
under  nucleation  of  only  a  leading  partial  dislocation.  Nucleation  of  the  trailing 
partial  dislocation  is  then  impossible  since  the  theoretical  shear  strength  r)|‘CP  is 
reached  before  the  critical  nucleation  load.  Moreover,  below  £  <  prtt  <  £Wb  any 
dislocation  mechanism  disappears  and  failure  occurs  bv  homogeneous  shear.  The 
competing  mechanisms  are  illustrated  in  Figure  3-12. 

Numerical  estimates  of  the  length  scales  for  some  metals  based  on  ab-initio  data 
[126]  are  summarized  in  Table  3.2.  illustrating  that  £]:Hf  and  p2rit  are  0  (inn). 
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Figure  3-11:  Subplot  (a)  shows  generalized  stacking  fault  curve  (Morse  potential,  compar¬ 
ison  with  EAM  copper).  Wo  consider  3  atomic  (111)  layers  in  the  weak  plane  on  cither  side 
of  the  crack  plane.  Subplot  (b)  depicts  the  theoretical  shear  strength  against  (111)  plane 
shearing  in  the  direction,  for  FCC  and  HCP  stacking  (Morse  potential).  The  HOP  stacking 
is  treated  as  a  completely  faulted  FCC  structure,  representing  a  close  approximation  to  the 
local  (111)  planar  arrangement  after  nucleation  of  the  first  partial  dislocation. 


71 


ss 


- - - - - ► 

"  Slab  size  c  (A) 


Figure  3-12:  Illustration  of  the  three  regimes  of  pure  homogeneous  shear  (ss),  partial 
dislocation  nucleation  and  shear  (ps),  and  1st  partial  and  2ud  partial  dislocation  nucleation 
(pp)  events  on  the  slip  plane.  The  two  curves  in  the  figure  show  nucleation  stress  as  a 
function  of  £,  for  the  1st  and  2nd  partials.  The  horizontal  lines  show  the  theoretical  shear 
strengths  for  FCC  and  HCP  stacking.  Intersection  of  the  1st  partial  nucleation  curve  and 
the  theoretical  shear  strength  for  FCC  crystal  yields  and  intersection  of  the  2nd  partial 
nucleation  curve  and  the  theoretical  shear  strength  for  the  HCP  crystal  yields  £%rtt  . 


Material 

(in)  [n^i 

(111)1110] 

nr*  (A) 

&rU  (A) 

zr*  (A) 

er*  (A) 

Copper 

25 

21 

23 

51 

Aluminum 

13 

11 

12 

13 

Table  3.2:  Summary  of  the  two  critical  length  scales,  £flt  and  Tlt  for  Cu  and  Al,  con¬ 
sidering  two  different  crack/slip  systems,  both  with  coincident  crack  and  slip  planes.  The 
estimates  are  based  on  ab-initio  data  for  stacking  fault  energies  and  shear  strengths  [126], 
for  which,  is  not  available.  Here,  it  has  been  assumed  that,  t$cp  =  rftcc 
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3.8  Atomistic  simulations  of  the  size  confinement 
effect 

The  above  derivation  predicts  a  change  in  deformation  mechanism  at  crack  tips  as 
the  system  size  is  reduced  to  the  order  of  a  few  nanometers.  Here  we  report  a  series  of 
atomistic  molecular  dynamics  simulations  to  directly  show  this  effect.  Simulating  this 
behavior  in  a  molecular  dynamics  simulation  provides  us  with  the  atomistic  details 
of  change  in  mechanism  as  the  system  size  is  gradually  changed. 

We  develop  a  simple  atomistic  system  with  coupled  interatomic  interactions  de¬ 
fined  by  Morse  [54]  and  harmonic  potentials  (Figure  3-10).  We  consider  a  two- 
dimensional  system  with  a  crack  along  a  weak  interface  joining  two  harmonic  crystals. 
All  atomic  interactions  are  harmonic  except  for  a  thin  layer  on  both  sides  of  the  crack 
plane.  Here,  pairs  of  atoms  which  are  separated  by  the  plane  of  the  crack,  ahead  of 
the  crack  tip,  are  allowed  to  interact  through  a  simple  Morse  potential  which  allows 
bonds  to  break  and  shearing  to  take  place.  The  crack  region  and  crack  faces  arc 
defined  by  zero  interactions  across  the  crack  plane  from  the  edge  of  the  sample  to  the 
crack  tip.  Harmonic  interactions  define  the  interactions  between  pairs  of  atoms  on 
the  same  side  of  the  shear  plane  on  both  sides  of  the  crystal.  The  choice  of  these  sim¬ 
plistic  potentials  is  dictated  by  our  goal  to  focus  on  fundamental  concepts  of  metallic 
systems,  that  is,  their  ability  to  undergo  localized  or  homogeneous  shear  deformation 
along  specific  slip  planes  for  different,  crystal  sizes.  The  use  of  harmonic  potential 
interactions  between  atoms  on  the  same  side  of  the  crack  plane,  allows  us  to  ensure 
existence  of  a  single  slip  plane  to  ensure  that  twinning  is  impossible  [127].  This  atom¬ 
istic  model  resembles  the  case  considered  in  the  theoretical  calculations.  Even  though 
wc  have  chosen  a  specific  thin  strip  geometry,  our  model  includes  the  key  character¬ 
istics  of  the  conditions  that  define  the  onset  of  plasticity,  that  is,  the  competition 
between  release  of  elastic  energy  and  the  energy  required  to  nucleate  partial  disloca¬ 
tions.  Since  plasticity  can  in  general  be  described  based  on  this  framework,  our  model 
provides  generally  valid  results,  for  cracks  of  any  size  under  geometric  confinement 
regardless  of  the  details  of  boundary  conditions. 
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We  emphasize  that  this  model  system  is  deliberately  designed  to  be  such  a  simple 
representation,  and  that  its  simplicity  is  not  a  due  to  the  lack  of  the  ability  to  build 
a  more  complex  representation  of  nanoscale  plasticity.  Further,  our  atomistic  model 
is  not  designed  to  make  quantitative  predictions.  Instead,  our  model  is  developed  to 
capture  the  most  significant  physical  quantities  and  processes  involved  in  dislocation 
nuclcation  -  clastic  energy  stored  in  the  bulk  (harmonic  potential)  released  to  induce 
localized  shear  in  a  slip  plane  (Morse  potential).  Thus,  for  the  specific  case  considered 
here,  our  atomistic  model  does  not  have  limitations  compared  with  EAM  potentials 
[56].  Similar  concepts  of  mixed  pair  potentials  (harmonic  potentials  of  different  types) 
have  been  used  earlier  to  study  hyperelasticity  at  crack  tips  in  dynamical  fracture 
[108], 

The  Morse  and  harmonic  potentials  used  in  the  atomistic  model  arc  defined  in 
section  2.1.2  and  reproduced  here: 

vab°rse  O’)  =  £ab  [{1  -  exp  (~aab  (r  -  cxq6))}2  -  l]  .  (3.8) 

and, 

VJT" M  =  \k.t  (r  -  r()„bf .  (3.9) 

The  variables  V*blorsea,nd  V^arm  are  the  potential  energies  for  a  and  b  atom  type 
interactions  of  Morse  and  harmonic  character  respectively.  eab.  aab,  aab  are  Morse 
parameters  and  kab,  r0ab  are  harmonic  potential  parameters  for  a-b  atom  type  inter¬ 
actions.  We  choose  sab=  0.2926  eV,  aab=  1.7866  A-1,  trQ&=2.7110  A  with  rcut=5.1  A 
for  the  Morse  potential;  and  kab=  20.0  eV/A,  r()ab=  2.675  A  (only  nearest  neighbor 
harmonic  interactions).  The  Morse  potential  parameters  lead  to  yus  =  200  mJ/in2 
and  7 Sy  =  40  m.J/m2  (values  for  Cu),  with  rth=4.25  GPa  (FCC)  and  3.25  GPa  (HCP) 
(Figure  3-ll(a)).  The  spring  constant  k  of  the  harmonic  potential  controls  the  clastic 
strain  energy  stored  in  the  bulk  crystal  under  shear  strain  and  is  set  to  a  value  such 
that  the  harmonic  bulk  modulus  is  ten  times  the  size  of  the  pure  Morse-potential  bulk 
modulus. 

The  system  considered  here  has  a  length  of  3.000  A  and  a  thickness  of  30  A, 
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whereas  the  height  parameter  £  is  varied  between  30  A  and  60  A  to  investigate  size 
effects.  The  crack  length  a  =  300  A  in  all  simulations.  The  system  is  loaded  in 
shear  by  displacement  boundary  conditions  at  the  top  anti  bottom.  A  niicrocanonical 
ensemble  with  a  time  step  of  2.5  fs  and  a  low  starting  temperature  at  0.1  K  is  used. 
The  shear  strain  rate  is  3.5E7  sec-1  (for  all  systems). 

3.9  Results  and  Discussion 

We  carry  out  a  systematic  study  of  the  deformation  mechanisms  while  varying  £. 
First,  we  study  the  shear  activity  at  the  slip  plane  under  the  lateral  shear  loading 
using  an  atomic  slip  vector  analysis  [91].  The  change  in  the  ^-component  of  the  slip 
vector  (measured  500  A  ahead  of  the  crack  tip)  versus  the  simulation  time  is  shown  in 
Figure  3-13,  measuring  the  relative  shear  displacement  of  the  upper  and  lower  part. 
We  observe  two  distinct  shear  events  for  each  slab  size.  Notably,  the  character  of 
the  shear  events  changes  from  a  smooth  continuous  profile  at  smaller  slab  sizes  to  a 
sudden  jump  at  larger  dimensions. 


Figure  3-13:  Variation  of  x-component  of  slip  vector  in  the  x-direction  with  bulk  shear 
strain,  at  a  distance  of  500  A  ahead  of  the  crack  tip.  The  curves  show  slip  vector  plots  for 
various  slab  heights. 


To  elucidate  the  atomistic  details  of  the  deformation  mechanism,  we  plot  the 


variation  of  the  .^-component  of  the  slip  vector  across  the  entire  slip  plane  in  the  re¬ 
direction  at  times  close  to  the  particular  shear  events  (Figure  3-14  for  the  first.  Figure 
3-15  for  the  second).  Comparing  these  plots  with  the  results  shown  in  Figure  3-13 
reveals  that  the  jagged  jumps  indeed  represent  the  emergence  of  partial  dislocations 
from  the  crack  tip,  and  the  smooth  transitions  represent  uniform,  homogeneous  shear 
across  the  whole  slip  plane. 

This  agrees  with  the  notion  that  each  shear  event  is  a  competition  between  nu- 
cleation  of  a  dislocation  at  the  crack  tip  and  a  homogeneous  shear  event  that  is  not 
localized  at  the  crack  tip.  After  the  first  partial  or  uniform  shear  event,  the  slip  plane 
has  a  stacking  fault  region  and  then  the  second  partial/uniform  shear  event  has  mod¬ 
ified  energy  barriers  to  overcome  because  of  the  change  in  stacking  sequence  ahead 
of  the  crack  tip.  If  the  second  event  is  a  non-local  homogeneous  shear  event  on  the 
slip  plane,  the  modified  barrier  is  the  energy  to  shear  planes  in  a  HCP  configuration; 
whereas  if  it  is  nucleation  of  the  second  partial,  the  modified  barrier  can  be  treated  as 
a  combination  of  a  reduced  energy  barrier  due  to  presence  of  a  higher  energy  stacking 
fault  and  also  the  stress  shielding  effect  of  the  presence  of  the  first  partial  in  the 
vicinity  of  the  crack  tip. 

The  results  shown  in  Figures  3-14  and  3-15  confirm  the  predictions  of  the  the¬ 
oretical  analysis.  For  the  smallest  slab  size  (37  A),  the  two  shear  events  represent 
homogeneous  shear  along  the  [112]  directions  on  the  slip  plane,  as  is  clearly  visible  in 
Figure  3-14(a)  and  3-15(a).  As  the  size  is  increased,  the  first  slip  event  changes  into 
nucleation  of  a  partial  dislocation  at  the  tip.  as  is  shown  in  Figure  3-14(b).  However, 
the  second  shear  event  remains  nonlocalized,  revealing  formation  of  a  shear  displace¬ 
ment  far  ahead  of  the  tip  (Figure  3-15(b)).  At  larger  slab  heights,  both  first  and 
second  shear  events  are  identified  as  partial  dislocations  nucleating  at  the  crack  tip. 
forming  a  complete  dislocation  (Figure  3-14(c)  and  3-15(c)). 

We  note  that  the  critical  length  scales  associated  with  the  nucleation  of  the  first 
and  second  partials  are  very  similar  for  the  (111)[112]  system  (Table  3.2).  Thus  the 
regime  of  a  single  partial  dislocation  (:ps;  in  Figure  3-12)  is  extremely  difficult  to 
observe.  This  is  reflected  by  the  results  shown  in  Figure  3- 15(b)  that,  while  providing 
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Figure  3-14:  Variation  of  x-componcnt  of  slip  vector  along  the  slip  plane  x-dircction  from 
the  crack  tip  to  the  slab  end  (first  shear  event),  for  slab  sizes  (37  A.  50  A  and  57  A).  The 
numbers  in  the  legend  show  the  corresponding  applied  strain  values.  The  results  show  a 
transition  from  uniform  shear  across  the  slab  length  ((a))  to  partial  dislocation  nucleation 
at  crack  tip  ((b)  and  (c)). 
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Figure  3-15:  Variation  of  the  x-component  of  the  slip  vector  along  the  slip  plane,  (a) 
reveals  uniform  shearing  across  the  slab  length,  (b)  reveals  that  shear  deformation  occurs 
also  ahead  ol  the  crack  tip  (marked  by  the  circle),  indicating  that  the  deformation  mecha¬ 
nism  is  not  controlled  by  the  crack  tip  stress  field  but  occurs  homogeneously  along  the  weak 
plane,  (c)  shows  nucleation  of  the  trailing  partial  dislocation. 


78 


some  evidence  for  nonlocal  slip,  do  not  clearly  reveal  homogeneous  shear  as  observed 
in  Figures  3-14(a)  and  3-15(a). 

In  summary,  our  theoretical  and  numerical  analysis  revealed  that  there  exist  fun¬ 
damental  length  scales  that  depend  only  on  material  parameters,  controlling  the  de¬ 
formation  mechanism.  The  theoretical  calculations  show  us  that  these  characteristic 
length  scales  are  on  the  order  of  a  few  nanometers  for  metals  such  as  Cu  and  Al. 
Moreover,  theoretical  calculations  for  inclined  slip  planes  at  the  crack  tip,  that  is, 
not  lying  in  the  crack  plane,  show  a  characteristic  length  scale  for  dislocation  nucle- 
ation  that  would  be  larger  by  a  factor  of  «2  than  this  simple  case  [111],  bringing 
the  numerical  estimates  to  the  range  of  5  to  10  urn  (based  on  the  estimates  shown 
in  Table  3.2).  Considerations  of  crystal  anisotropy  change  the  characteristic  length 
scale  further  by  ±  50%  [128].  These  change  the  characteristic  length  scales  while 
maintaining  the  same  order  of  magnitude. 

A  crack  tip  stress  held  rises  from  a  discontinuity  in  interactions  along  the  crack 
tip  line  across  a  plane,  giving  rise  to  a  stress  field  ~  1/y/r  for  a  semi-infinite 
crack  under  far-field  uniform  K\\  loading.  Whereas  our  study  has  been  limited  to  a 
cracked  thin  strip  geometry,  one  can  extend  these  models  to  other  cases,  e.g.  to  rigid- 
compliant  interfaces  in  composites  or  grain  boundaries  in  nanocrystalline  materials. 

Our  analysis  [129]  is  directly  applicable  in  the  design  of  bone-inspired  metal¬ 
lic  nanocomposites,  covered  earlier.  Figure  3-1G  illustrates  the  application  of  the 
size  effect  discussed  in  this  section  for  the  analysis  of  the  flow  stress  of  a  bioin¬ 
spired  nanocomposite  discussed  earlier  [106].  Figure  3-16(a)  shows  the  geometry  of 
a  nanocomposite  composed  of  a  soft,  ductile  matrix  and  hard,  brittle  platelets.  The 
application  of  uniaxial  tensile  load  reveals  that  large  shear  stresses  are  transferred 
between  the  platelets  as  shown  in  Figure  3-16(b),  leading  to  nucleation  of  disloca¬ 
tions  in  the  matrix  phase  during  plastic  deformation.  Therefore  the  characteristic 
dimension  £  (characterizing  the  thickness  of  the  ductile  matrix  material)  controls  the 
flow  stress  level,  and  at  small  dimensions  below  the  characteristic  length  scales  it 
is  expected  that  dislocation  based  deformation  breaks  down.  Figure  3- 16(c)  depicts 
an  analysis  of  the  flow  stress  as  a  function  of  the  building  block  size,  for  a  different 
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Figure  3-16:  Illustration  of  the  application  of  the  size  effect  discussed  here  for  the  analysis 
of  the  flow  stress  of  a  bioinspired  nanocomposite  [106].  Subplot  (a)  shows  the  geometry 
of  a  nanocomposite,  also  shown  in  Figure  3-1.  The  application  of  uniaxial  tensile  load 
reveals  that  large  shear  stresses  are  transferred  between  the  platelets  (subplot  (b)),  leading 
to  nucleation  of  dislocations  in  the  matrix  phase  during  plastic  deformation.  Therefore  the 
characteristic  dimension  f  controls  the  flow  stress  level,  and  at  small  dimensions  below  the 
characteristic  length  scales  it  is  expected  that  dislocation  based  deformation  breaks  down. 
Subplot  (c)  depicts  an  analysis  of  the  flow  stress  as  a  function  of  the  building  block  size, 
for  a  Ni-Al  nanocomposite  [105]  possessing  the  same  design  morphology,  showing  that  the 
flow  stress  increases  with  a  reduction  of  the  building  block  size  (notably,  the  strengthening 
behavior  agrees  well  with  the  predicted  scaling)  then  reaches  a  maximum  at  a  critical 
building  block  size,  followed  by  a  decay.  This  transition  corresponds  to  the  breakdown  of 
dislocation  activity  in  favor  of  intcrfacial  slip  (between  Ni  and  A1  particles),  as  is  revealed 
by  the  analysis  of  interfacial  slip  (green  curve  in  subplot  (c),  transition  from  regime  I  to 
regime  II). 
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Ni/Al  nanocomposite  with  the  same  design  morphology  [105],  showing  that  indeed 
the  flow  stress,  increases  with  reduction  of  £,  then  reaches  a  maximum  at  a  critical 
building  block  size,  followed  by  decay.  This  transition  corresponds  to  the  breakdown 
of  dislocation  activity,  as  is  revealed  by  the  analysis  of  intcrfacial  slip  (green  curve, 
transition  from  regime  I  to  regime  II).  Notably,  the  strengthening  behavior  agrees  well 
with  the  predicted  scaling  rfrit  ~  1  f  \f^,  providing  further  evidence  that  the  simple 
model  captures  the  essential  features  of  the  flow  stress  hardening  process. 

In  nanocrystals,  the  grain  boundaries  are  sites  of  weaker  stress  concentration  than 
crack  tips,  and  there  exist  several  other  accommodating  mechanisms.  Our  system  was 
deliberately  designed  so  that  it  only  contains  two  failure  mechanisms;  realistic  systems 
contain  other  failure  mechanisms  (grain  boundary  sliding,  twinning,  diffusion).  Thus, 
in  nanocrystals,  change  in  plasticity  mechanism  due  to  geometric  confinement  can  be 
considered  to  be  only  one  of  several  length  scales  of  importance  in  the  competition 
between  different  deformation  mechanisms  [123,  130,  131,  132], 

Since  these  other  geometries  feature  a  weaker  stress  concentration,  our  analysis 
provides  a  lower  bound  for  the  critical  length  scales.  This  can  be  explained  by 
considering  the  scaling  relation  of  the  change  of  the  characteristic  length  scales.  The 
condition  for  nuelcation  of  a  dislocation  can  be  generally  written  as  GfIiti  =  G.  Since 
G  ~  Kjj  ~  £.  a  change  in  Kn  to  ctKn  leads  to  a  change  of  £  ~  Gfj\/a2.  Therefore 
for  a  <  1  the  characteristic  length  scales  £?"*  increase.  Thus,  while  the  basic  physical 
mechanisms  are  captured  appropriately  in  our  simple  model,  for  many  applications 
the  critical  length  scales  may  be  several  times  larger  than  those  reported  here. 

Our  theoretical  system  was  deliberately  designed  so  that  it  only  contains  two  fail¬ 
ure  mechanisms.  Realistic  systems  contain  other  failure  mechanisms  (grain  boundary 
sliding,  twinning,  diffusion).  This  may  explain  why  transitions  of  mechanisms  in 
experiments  and  simulations  of  nanocrystalline  metals  have  been  observed  at  larger 
scales  than  those  listed  in  Table  3.2.  MD  simulations  [105,  130,  131,  133]  and  ex¬ 
periments  [123]  of  nanocrvstalline  metals  have  shown  that  that  at  grain  sizes  below 
~20  nm,  all  dislocation  plasticity  arises  at  grain  boundaries.  Further,  this  gives  rise 
to  twins  and  stacking  faults  across  the  grains,  that  is,  complete  dislocations  are  not 
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observed  at  these  grain  sizes  [130,  131,  132.  133]. 

The  partial  dislocations  that  give  rise  to  these  features  are  emitted  from  grain 
boundary  dislocation  networks  close  to  triple  junctions  [130] .  These  are  possible  loca¬ 
tions  where  high  energy  and  low  energy  grain  boundaries  meet,  and  stress  relaxations 
along  particular  grain  boundaries  leads  to  stress  concentration  at  the  junctions.  The 
atom  packing  and  thus  interatomic  bonding  is  quite  different  between  these  sections 
of  the  grain  boundaries  which  produce  a  stress  concentration  field  at  the  intersection 
of  these  regions  in  the  plane  of  the  grain  boundary  when  a  shear  stress  is  applied 
across  the  two  grains  sharing  the  boundary.  The  change  in  bonding  arrangement 
is  not  as  sharp  as  for  a  crack  tip.  Thus  the  stress  field  falls  off  much  slower  as  a 
function  of  distance  to  the  grain  boundary  than  it  is  the  case  for  a  sharp  crack.  This 
implies  that  our  calculations  are  a  lower  bound  on  the  characteristic  crystal  sizes 
that  have  change  in  mechanism.  In  agreement  with  this  notion,  at  sizes  below  ~10 
nm  both  experiments  and  simulations  show  breakdown  of  any  dislocation  activity 
either  from  grain  boundaries  or  the  bulk.  All  plastic  deformation  is  accounted  for  by 
grain  boundary  sliding,  atomic  shuffling  or  grain  boundary  diffusion  as  accommodat¬ 
ing  mechanisms  [130,  131].  This  shows  a  fundamental  application  of  our  length  scales 
which  limit  dislocation  nucleation  to  the  problem  of  change  in  plasticity  mechanism 
as  grain  size  is  reduced  in  the  nanometer  scale.  Even  though  the  model  reported 
here  does  not  provide  a  quantitative  link,  it  leads  to  an  explanation  of  the  underlying 
physical  mechanisms. 

Similar  considerations  apply  for  the  mechanics  of  thin  film  systems.  It  has  been  re¬ 
ported  that  dislocation  activity  is  absent  in  nanocrystalline  thin  films  of  average  grain 
size  of  10  20  nm  grain  size  under  tensile  deformation  [120,  133].  Earlier  MD  studies 
have  also  shown  that  plasticity  from  threading  and  parallel  glide  plane  dislocations 
is  seen  to  be  absent  at  a  length  scale  below  25  nm  [133]. 

Our  theoretical  analysis  has  also  used  in  the  design  of  bone-inspired  metallic 
nanocomposites.  Figure  3-16  illustrates  an  application  of  the  size  effect  discussed 
in  this  study  for  the  analysis  of  the  flow  stress  of  a  bioinspired  nanocomposite.  Fig¬ 
ure  3-lG(a)  shows  the  geometry  of  a  nanocomposite  composed  of  a  soft,  ductile  matrix 
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and  hard,  brittle  platelets.  The  application  of  uniaxial  tensile  load  reveals  that  large 
shear  stresses  are  transferred  between  the  platelets  as  shown  in  Figure  3-16(b),  leading 
to  nucleation  of  dislocations  in  the  matrix  phase  during  plastic  deformation.  There¬ 
fore  the  characteristic  dimension  £  (characterizing  the  thickness  of  the  ductile  matrix 
material)  controls  the  flow  stress  level,  and  at  small  dimensions  below  the  charac¬ 
teristic  length  scales  it  is  expected  that  dislocation  based  deformation  breaks  down. 
Figure  3-16(c)  depicts  an  analysis  of  the  flow  stress  as  a  function  of  the  building 
block  size,  showing  that  indeed  the  flow  stress,  increases  with  reduction  of  £,  then 
reaches  a  maximum  at  a  critical  building  block  size,  followed  by  decay.  This  transition 
corresponds  to  the  breakdown  of  dislocation  activity,  as  is  revealed  by  the  analysis 
of  interface  slip  (green  curve,  transition  from  regime  I  to  regime  II).  Notably,  the 
strengthening  behavior  agrees  well  with  the  predicted  scaling  t[''iL  ~  lj providing 
further  evidence  that  the  simple  model  captures  the  essential  features  of  the  flow 
stress  hardening  process. 

To  the  best  of  our  knowledge,  this  study  is  the  first  to  systematically  show  intrin¬ 
sic  length  scale  limitations  to  plasticity,  along  with  direct  confirmation  by  atomistic 
simulation.  A  related  behavior  has  been  observed  in  simulations  of  failure  of  thin 
metal  nanowires  (diameter  f»15  A)  under  tension,  where  a  competition  between  ho¬ 
mogeneous  shear  and  dislocation  nucleation  is  controlled  by  the  energy  of  slip  across 
cross-section  of  nanowire  and  the  self-energy  of  the  dislocation  [134].  The  analysis 
reported  here  is  similar  to  earlier  work  in  brittle  materials  [33],  where  a  single  critical 
length  scale  has  been  proposed  below  which  no  brittle  fracture  mechanism  can  oc¬ 
cur  (the  ‘flaw-tolerance’  length  scale).  Interestingly,  for  plasticity,  two  critical  length 
scales  exist,  provided  that  full  dislocations  are  split  up  into  two  partial  dislocations 
as  it  occurs  in  many  FCC  metals. 


3.10  Conclusions 

Our  analysis  reveals  that  there  exist  fundamental,  intrinsic  length  scales  that  depend 
only  on  material  parameters  and  the  particular  geometry  that  control  the  plastic  de- 
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formation  mechanism  and  strength  properties  in  metallic  nanocomposites.  Our  study 
may  be  vital  in  the  analysis  and  design  of  new  nanomaterials,  for  instance  bioinspired 
nanocomposites,  providing  concrete  design  suggestions  and  associated  mechanistic 
models.  Specifically,  in  agreement  with  earlier  findings  [135,  33]  we  find  that  the  use 
of  elongated  platelets  of  high  aspect  ratio  and  staggered  arrangement  of  platelets  for 
optimal  load  transfer,  and  control  of  spacing  between  layers  of  platelets  arc  critical 
factors  in  strengthening  the  material.  We  also  observe  that  the  strength  of  the  matrix- 
reinforcement  interface  determines  the  optimal  size  of  the  second  phase  platelets  at 
which  the  maximal  flow  strength  is  observed. 

To  the  best  of  our  knowledge,  unlike  earlier  studies  of  similar  geometries,  here  wc 
have  for  the  first  time  explored  a  much  larger  parameter  space  and  provided  a  detailed 
investigation  of  associated  mechanical  properties  and  their  relationship  to  the  under¬ 
lying  nanostructural  design.  This  provides  details  insight  into  structure-property 
relationships  in  this  class  of  materials  from  a  bottom-up  atomistic  perspective  that 
explicitly  includes  dislocation  mechanics  and  interfacial  phenomena. 

We  summarize  the  main  results  of  the  study  of  dependence  of  strength  of  bone- 
inspired  metallic  nanocomposites  on  their  design: 

1.  Reinforcing  phase  platelets  arranged  in  staggered  rows  with  the  stagger  perpen¬ 
dicular  to  the  direction  of  loading  lead  to  optimal  load  transfer  in  the  nanocom¬ 
posite.  A  stagger  of  0.25  to  0.75  leads  to  a  characteristic  tension-shear  load 
transfer  in  platelets /matrix  that  maximizes  the  flow  stress. 

2.  Reinforcing  phase  platelets  with  large  aspect  ratio  with  the  long  edge  in  the 
direction  of  loading  lead  to  higher  flow  stress.  Beyond  a  critical  aspect  ratio 
(larger  than  4  for  materials  considered  here),  the  flow  stress  saturates,  and  this 
parameter  has  no  further  effect. 

3.  The  spacing  between  platelets  perpendicular  to  the  direction  of  loading  (trans¬ 
verse  spacing)  has  a  Hall-Petch  like  effect  on  flow  stress  with  the  spacing  mim¬ 
icking  the  grain  diameter. 
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4.  The  volume  fraction  of  the  Second  phase,  while  keeping  the  aspect  ratio  and 
platelet  spacing  constant,  does  not  have  a  significant  effect  on  flow  stress.  This 
result  provides  strong  evidence  for  the  significance  of  the  specific  geometry,  and 
suggests  that  nanoscale  control  of  structure  and  shape  of  platelets  is  crucial  to 
maximize  material  performance. 

5.  A  comparison  of  the  effect  of  various  geometrical  effects  on  flow  stress  shows 
that  a  factor  of  1.94  improvement  over  the  unreinforced  matrix  flow  strength  is 
possible  by  systematically  engineering  the  nanostructured  design. 

G.  The  interfacial  strength  between  matrix  and  platelets  is  a  key  parameter  that 
determines  the  importance  of  sliding  and  decohesion  as  deformation  mechanisms 
at  small  spacing  and  platelet  sizes.  This  leads  to  a  peak  flow  strength  for  the 
nanocomposite  as  a  function  of  platelet  size  and  spacing. 

Our  analysis  also  reveals  that  there  exist  fundamental,  intrinsic  length  scales  that 
depend  only  on  material  parameters  and  the  particular  geometry  that  control  the 
plastic  deformation  mechanism  in  small  crystals  under  confined  conditions.  These 
characteristic  length  scales  separate  regimes  of  no  dislocation  activity,  partial  dislo¬ 
cation  plasticity,  and  complete  dislocation  plasticity  at  a  crack  tip  in  ductile  metals. 
We  have  confirmed  this  effect  by  direct  atomistic  simulation  of  a  model  system  under 
shear  mode  II  loading  with  coincident  crack  and  slip  planes.  Wc  have  also  shown 
the  applicability  of  this  concept  to  the  interfacial  effect  on  strength  wc  observe  in 
the  metallic  nanocompositcs.  Confined  ductile  phases,  in  such  materials,  will  show  a 
transition  to  homogeneous  shear  based  plasticity  below  a  critical  length  scale.  This 
could  provide  important  guidance  for  the  optimal  design  of  such  nanocomposites,  as 
the  ductile  phase  will  fail  at  their  theoretical  strength  and  any  further  reduction  in 
the  critical  dimension  will  not  increase  the  failure  strength. 

The  study  reported  here  illustrates  that  universally  observed  biological  nanos¬ 
tructures  might  provide  effective  design  solutions  for  structurally  strong  materials. 
Thereby,  the  design  of  a  material’s  nanostructure  is  a  crucial  element  in  achieving 
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superior  mechanical  properties.  By  utilizing  functional  building  blocks  such  as  met¬ 
als,  ceramics  or  polymers  it  might  be  possible  to  add  additional  functionality  to  the 
materials,  beyond  structural  mechanical  properties.  By  the  control  of  interfacial  prop¬ 
erties,  e.g.  by  sensitive  polymers  that  respond  to  external  cues  such  as  magnetic  fields, 
light  or  chemical  environment  [136],  it  might  be  possible  to  engineer  novel  responsive, 
active  and  tunable  materials  with  substantial  variations  in  mechanical  properties. 

The  use  of  metals  in  the  design  of  bioinspired  structures  is,  however,  a  costly 
proposition,  due  to  the  economic  cost  of  the  raw  materials.  The  design  of  bioinspired 
structures  would  be  economically  much  more  feasible,  if  it  could  be  undertaken  with 
very  cheap  raw  materials.  Silica,  found  abundantly  in  sand,  is  one  such  cheap,  however 
structurally  'poor’  material.  Is  it  possible  to  use  silica  as  the  building  block  of  a 
bioinspired  structure,  which  will  provide  excellent  mechanical  properties?  In  the  next 
chapter,  wc  turn  for  inspiration  to  the  exoskeleton  of  diatoms,  a  silica-based  biological 
material.  Wc  present  the  design  and  mechanical  properties  of  diatom-inspired  silica 
nanostructures. 
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Chapter  4 


Ductility  enhancement  through 
diatom-inspired  nanoporous  silica 
design 


Silica  is  one  the  most  abundant  minerals  in  the  earth's  crust,  known  for  its  hardness 
and  brittle  fracture  behavior.  Silica  is  its  amorphous  glassy  form  or  crystalline  quartz 
form  is  considered  a  prototype  of  a  perfectly  brittle  ceramic  material,  with  little  to 
no  plastic  deformation  prior  to  fracture,  and  thus  following  the  Griffith’s  criterion  of 
fracture  [137,  138].  In  biology,  silica  structures  have  been  observed  to  be  assembled  up 
from  the  nanoscale  in  honeycomb  and  porous  form  [139].  For  example,  nature  shows 
the  use  of  porous  silica  structures  in  the  exoskeleton  of  diatoms  [140,  141,  142.  143]. 
In  diatoms,  a  porous  hierarchical  structure  dominates  the  landscape  of  the  cell  wall 
and  encompasses  intricate  patterns  that  are  highly  varied  and  ordered,  as  shown  in 
Figure  1-2,  and  reach  all  the  way  down  to  the  nanoscale.  A  marine  diatom  species 
( Concinodicus  sp.),  seen  in  Figure  1-2,  possesses  a  silica-based  cxoskclcton  (called 
frustulc-cxtcrnal  surface  of  the  diatom  at  a  micron  scale)  made  up  of  porous  parts 
arranged  in  a  hierarchical  fashion,  with  areola  pores,  the  internal  surface  of  the  diatom 
(at  micron  length  scale);  the  2nd  central  porous  layer,  the  cribrum  (with  pores  at  a 
sub-micron  length  scale);  and  the  cribellum,  the  external  porous  layer  (with  pores 
at  a  nanometer  length  scale).  Beyond  the  key  biological  functions  of  the  cell  wall 
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lie  important  mechanical  functions,  such  as  preventing  virus  penetration,  protecting 
diatoms  from  the  jaws  of  predators,  or  even  protection  from  digestion  in  some  cases 
[144.  145].  It  is  intriguing  that  a  material  that  is  so  brittle  in  its  engineering  form,  is 
used  by  nature  as  the  major  constituent  of  the  protective  casing  of  these  species. 

The  question  of  the  mechanical  properties  of  these  silica  structures  is  key  to  un¬ 
derstanding  their  use  in  the  exoskeletons  [39].  Some  recent  experimental  studies  that 
revealed  the  mechanical  properties  of  diatom  shells  were  covered  in  Chapter  1,  section 
secfiiierarchyeffectfailure.  Briefly,  they  have  been  found  to  possess  high  strength,  and 
carry  large  reversible  elastic  strains.  Also  large  variations  in  mechanical  properties 
have  been  observed  to  be  influenced  by  pore  size,  pore  distance,  porosity  and  under 
different  biomineralization  processes  [40,  41]. 

In  accordance  with  our  aim  in  Chapter  1,  we  thus  want  to  study  the  effect  of 
(a)  nanostructuring,  and,  (b)  the  use  of  hierarchical  assembly,  on  silica  material 
design.  We  study  in  this  chapter  firstly,  the  effect  of  nanostructuring  on  silica,  and 
the  mechanical  properties  of  nanoporous  diatom- inspired  silica  structures. 


4.1  Background  on  nanoscale  silica  structures 

Nanoporous  silica  structures  have  been  also  designed  in  the  lab  [146,  147].  and  have 
great  utility  as  catalysts  and  adsorption  media.  They  could  also  be  used  as  tem¬ 
plates  for  assembling  other  materials  into  nanostructures.  However,  their  mechanical 
properties  have  not  been  investigated,  and  it  remains  unknown  how  they  perform 
mechanically,  in  particular  under  extreme  loading  or  deformation.  This  is  important 
for  the  reliability  of  the  application  devices,  where  issues  such  as  fracture  proper¬ 
ties  are  essential.  A  detailed  analysis  of  the  mechanics  of  nanoporous  silica  struc¬ 
tures  can  also  be  used  to  discover  design  criteria  to  maximize  their  load-carrying  and 
fracture-resistant  capacity,  and  open  up  new  applications  for  these  structures  in  novel 
materials. 

One-dimensional  (ID)  struct  ures  of  silica,  silica  nanowires,  have  also  been  under 
investigation  experimentally  and  theoretically  for  possessing  potential  applications 
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as  structural  building  blocks  in  nanoelectronics  and  optics.  Precise  measurement 
of  their  mechanical  properties  is  important,:*  as  mechanical  failure  of  these  blocks 
may  lead  to  failure  of  the  entire  devices.  Several  experiments  show  that  there  exists 
a  size-dependent  brittlc-to-ductile  transition  in  these  nanowires.  Amorphous  silica 
nanowircs  of  50-100  run  diameters,  obtained  using  chemical  vapor  deposition  tech¬ 
niques.  have  been  shown  to  exhibit  brittle  fracture  under  bending  loads  [148].  However 
silica  nanowires  of  diameters  down  to  20  nm,  created  using  physical  taper-drawing, 
show  pristine  smooth  surfaces  and  great  deform  ability,  with  the  ability  to  be  assem¬ 
bled  into  spiral  coils,  nanorings  and  nanoloops  [149],  Also,  crystalline  SiC.  another 
brittle  ceramic,  has  been  shown  to  exhibit  plastic  deformation  for  nanowires  with 
diameters  below  100  nm.  under  bending  loads  through  lattice  disordering  and  amor- 
phization.  Silicon  nanowires  with  diameters  below  60  nm  show  large  plastic  defor¬ 
mation  before  fracture,  through  emission  of  dislocations,  amorphization  and  necking 
[150). 


An  understanding  of  silica  behavior  at  the  nanoscale  is  thus  critical  to  under¬ 
standing  the  mechanics  of  nanoporous  silica  structures  and  their  role  in  the  structure 
of  diatom  exoskeletons,  and  also  the  mechanics  of  silica  nanowires.  Though  there 
exists  an  abundance  of  literature  on  their  optical  properties,  the  mechanical  proper¬ 
ties  of  silica  structures  at  the  nanoscale,  whether  they  are  nanowires  or  nanoporous, 
have  been  difficult,  to  measure  experimentally  due  to  problems  in  fabricating  struc¬ 
tures  with  no  surface  defects  (smooth  surfaces)  and  problems  in  loading  setup,  with 
some  exceptions  [149,  150].  Atomistic  simulations  with  accurate  descriptions  of  silica 
behavior  arc  a  critical  way  to  obtain  insight.  Previous  work  by  the  authors  oil  the 
modeling  of  the  behavior  of  nanoporous  silica  structures  using  atomistic  modeling 
has  shown  increased  ductility  with  a  reduction  in  the  constituent  silica  strut  size, 
and  a  change  of  fracture  behavior  from  brittle  crack  propagation  to  ductile  failure 
with  size  scaling  [151].  In  this  work,  we  use  atomistic  simulations  and  theoretical 
nano-mechanical  analysis  to  understand  the  enhanced  ductility  and  toughness  seen 
in  these  nanoporous  structures. 
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4.2  Design  parameters  of  the  nanoporous  silica  struc¬ 
tures 

Here  we  follow  the  structural  framework  of  diatoms  closely,  such  that  it  resembles 
their  porous  structure  at  the  smallest  hierarchy  level  of  the  cribellum,  with  pores  at  a 
nanometer  size  scale,  resulting  in  a  similar  model  system  of  silica  as  shown  in  Figure 
4-1.  The  nanoporous  structures  thus  designed  will  also  subsequently  be  referred  to  as 
nano-honeycomb  structures  due  to  their  2-dimensional  periodic  symmetry  and  simi¬ 
larity  with  macroscopic:  honeycombs.  The  widths  w  of  the  nano- honeycomb  structures 
arc  controlled  in  our  analysis,  and  range  from  5  A  to  72  A  (the  upper  size  is  limited 
by  computational  resources).  The  approach  pursued  here  is  guided  by  our  desire  to 
develop  a  general  model  system  in  which  we  can  test  the  effect  of  the  size  (geometric 
confinement)  of  the  nanostructure  on  the  bulk  material  behavior.  By  systematically 
varying  the  size  and  hierarchy  of  the  constituting  silica  nanostructure,  we  examine 
associated  mechanical  properties,  as  well  as  fracture  and  toughening  mechanisms, 
facilitated  through  a  series  of  molecular  dynamics  simulations. 


4.3  Materials  and  methods 

At  the  nanoscale  we  use  molecular  dynamics  simulations  to  study  the  mechanics  of 
the  nanoporous  silica  structures.  To  describe  the  inter-atomic  interactions,  we  use 
the  first-principles  derived  ReaxFF  force  field  [62]  (see  Chapter  1.  section  2.1.2). 
The  ReaxFF  description  for  Si-0  systems,  is  based  on  a  bond-length  bond-order 
description  and  fitted  to  density- functional  calculations  of  energy  landscapes  of  bond- 
distortion,  breaking  and  forming  events  of  various  Si-0  reactions.  A  variety  of  Si-0 
clusters  are  used  for  fitting  parameters,  as  also  energetics  of  bulk  crystalline  phases  of 
silicon  and  silica  under  tension  and  compression.  The  ReaxFF  potential  has  been  used 
successfully  in  predicting  fracture  phenomena  in  silicon  and  silica  [86,  104.  70,  152], 
and  interfacial  structure  at  silicon/silica  interfaces  [71.  153]. 
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n  =  72  A 


Figure  4-1:  Geometry  of  the  bioinspired  silica  structure,  and  setup  used  in  our  simulations. 
Panel  (a):  Three-dimensional  schematic  of  the  silica  nano-honeycomb  structure  (shown  on 
left),  with  periodic  boundaries  along  the  X,  Y.  and  Z  directions.  On  the  right,  the  geometry 
and  loading  of  the  individual  silica  struts  are  shown.  The  crystallographic  orientation  is  the 
same  for  all  silica  structures  considered.  The  arrows  indicate  tensile  load  applied  uniformly 
along  the  structure.  Panel  (b):  Initial  geometry  of  nano-honeycomb  structures  considered 
here,  illustrating  the  wall  width  ( w.  definition  indicated  in  one  of  the  structures)  variation 
in  the  geometry.  The  inset  shows  a  detailed  view  of  the  relaxed  surface  structure. 
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4.4  Deformation  of  nano-honeycomb  silica  struc¬ 


tures 


Here  we  present  our  analysis  on  the  effect  of  changing  the  wall  width  on  the  me¬ 
chanical  properties  of  the  silica  nano-lioneycomb(see  Figure  4- 1(a)  for  the  geometry 
considered).  As  shown  in  Figure  4-1.  the  honeycomb  wall  widths  are  varied  between 
w  =  5  A  to  72  A  for  the  structures. 

Figure  4-2  shows  the  stress-strain  response  of  the  nano-honeycomb  structure  for 
varying  widths  w.  The  structures  show  an  increase  of  deformation  in  the  plastic 
regime,  a  lower  modulus,  and  lower  maximum  stress  with  decreasing  wall  width. 
Even  though  silica  is  considered  a  brittle  material,  the  results  show  that  it  is  possible 
to  transform  it  into  a  ductile  system  for  small  nanoscale  wall  widths  which  reach  a 
maximum  failure  strain  of  ~120%. 


0  0.5  1  1.5 

Engineering  strain 


Figure  4-2:  Stress-strain  graphs  for  silica  nano-honeycomb  structure  for  all  sizes  (wall 
widths  range  from  w  =  5  A  to  72  A).  For  wall  widths  below  ~60  A,  we  observe  the  existence 
of  a  plastic  deformation  regime.  The  greatest  deformation  is  obtained  for  the  smallest  wall 
width  of  5  A.  Failure  mechanisms  are  characterized  by  crack  initiation  and  propagation,  or, 
shearing  followed  by  nanoscale  void  formation  and  coalescence. 


We  proceed  with  analysis  of  the  silica  meshes  as  shown  in  Figure  4-3,  which  shows 
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the  equivalent  von  Mises  stress  field  at  the  maximum  stress.  For  larger  wall  widths, 
high  stress  is  mostly  concentrated  on  the  surface  and  specifically  near  the  edges, 
thus  suggesting  possible  locations  for  crack  or  shear  nucleation.  However,  with  lower 
wall  widths,  the  stresses  become  relatively  homogeneous  throughout  the  structure. 


For  large  deformation,  the  void  shapes  gradually  change  from  a  rectangular  to  a 
hexagonal  one  for  decreasing  wall  widths,  and  can  be  clearly  seen  for  w  <  31  A. 


w  =  15  A 


w=  31  A 


Figure  4-3:  Von  Mises  stress  field  at  the  maximum  load,  for  different  nano-honeycomb 
wall  widths  (the  strain  value  at  which  the  snapshot  was  taken  is  indicated  in  the  plot).  For 
widths  smaller  than  s;31  A.  the  structure  at  the  maximum  load  becomes  hexagonal,  and 
the  stress  is  distributed  homogeneously  throughout  the  structure.  For  larger  wall  widths, 
high  stresses  are  concentrated  around  the  corners.  Moreover,  the  initial,  rectangular  shape 
of  the  structure  is  maintained.  In  order  to  improve  image  clarity,  we  only  show  the  stress 
values  associated  with  silicon  atoms  within  the  silica  system. 


The  analysis  sheds  some  light  on  the  remarkable  stress-strain  response  of  nano¬ 
honeycomb  structures  with  thin  wall  widths,  as  shown  in  Figure  4-2.  A  key  feature 
is  the  geometric  pattern  that  allows  large  deformations  to  be  accommodated  bv  the 
mesh  by  changing  from  a  rectangular  pattern  to  a  hexagonal  one  at  large  strains  (see, 
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e.g.  in  Figure  4-3).  specifically  for  wall  widths  below  31  A.  The  reason  for  these  very 
large  strains  without  failure  seems  to  be  due  to  the  more  homogeneous  distribution  of 
stresses  and  the  geometry  transformation  from  rectangular  to  a.  hexagonal  shape  for 
smaller  wall  widths.  In  the  next  few  sections,  we  develop  detailed  theoretical  models 
to  explain  and  predict  this  increase  in  ductility  with  decreasing  wall  size  precisely. 

4.5  Analysis  of  deformation  using  theoretical  mod¬ 
els 

The  atomistic  studies  show  the  effect  of  pore  size,  distribution  and  porosity  on  elastic 
modulus,  plasticity,  ductility  and  toughness  of  the  structures.  In  the  next  section, 
the  underlying  mechanisms  are  discovered  through  a  systematic  analysis  of  the  earlier 
simulation  results.  Also,  continuum  elasticity  theory  and  fracture  mechanics,  using 
atomistic  simulation  data,  are  used  to  make  theoretical  predictions  of  structural  stiff¬ 
ness  and  strength  for  varying  geometry.  The  theoretical  predictions  provide  design 
guidelines  for  how  to  make  nanoporous  materials  that  show  enhanced  ductility,  plastic 
flow,  and  toughness. 

4.6  Results  and  discussion 

We  study  the  mechanics  of  silica  nanoporous  structures  whose  morphology  is  shown 
in  Figure  4-4a.  The  structures  consist  of  rows  of  rectangular  pores  arranged  in  a 
staggered  fashion,  as  seen  in  several  natural  and  man-made  nanoporous  structures. 
The  structures  can  be  uniquely  identified  by  the  inter-pore  distance  and  the  size  of  the 
pores,  and  will  hereafter  be  referred  to  by  three  numbers,  (pore  distance  t,  pore  length 
Pi:  pore  width  pw)  (see  Figure  4-4a).  The  material  considered  is  silica  in  its  a-quartz 
polymorph,  which  is  the  stable  crystalline  form  at  room  temperature  and  pressure. 
The  structure  is  free  of  any  grain  boundaries  or  initial  defects.  All  structures  are 
periodic  in  all  3  directions,  and  loading  is  achieved  by  deformation  of  the  system  in 
the  Y  direction.  System  evolution  is  studied  through  molecular  dynamics  simulations 
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using  the  canonical  ensemble  with  the  Berendsen  thermostat  [47] .  A  time  step  of  0.2 
fs  is  used  to  advance  the  dynamics  and  a  deformation  strain  rate  of  1  x  1 0los“ 1  is 
applied.  All  simulations  are  carried  out  in  GRASP  [154].  a  parallelized  code  for  large 
scale  RoaxFF  simulations.  Figure  4-4(b)  shows  the  stress-strain  curves  for  some  of  the 
different  nano- honeycomb  structures  of  silica  that  were  analyzed.  There  is  a  marked 
difference  from  a  typical  macroscopic  brittle  ceramic  honeycomb  stress-strain  curve, 
which  shows  brittle  fracture  at  a  few  percent  strain  [139]. 


strain 


Figure  4-4:  Geometric  shape  classification  and  stress-strain  behavior  for  the  nano- 
honeycomb  silica  structures;  several  more  structures  showing  ductile  response  are  analyzed 
here,  apait  fiom  Figure  4-2.  (a)  shows  the  geometry  of  the  nano- honeycomb  structures 
under  consideration  with  3  independent  parameters  required  for  the  geometric  shape  clas¬ 
sification  (c.g.  t.  pi  and  pw),  (b)  show  stress  (oyy)-strain  curves  obtained  from  atomistic 
simulations  of  different  nano-honeycomb  structures.  The  legend  shows  the  classification  of 
the  structures,  which  is  shown  as  (t,  pi,  pw)  parameters  for  each  structure,  values  given  in 
A.  The  (323  A,  78  A.  30  A)  structure  shows  purely  brittle  fracture,  the  (56  A,  78  A.  30  A) 
structure  shows  simultaneous  plasticity  and  crack  propagation  in  different  struts,  the  rest 
of  the  structures  show  ductile  fracture. 


4.6.1  Elasticity 


The  nanostructure  can  be  considered  to  consist  of  a  collection  of  interlocking  struts, 
some  of  which  are  initial  horizontally  oriented  (length  in  the  X  direction)  and  the 
others  vertical  (length  in  the  Y  direction)  and  are  referred  to  as  such  in  the  rest  of 
this  study.  Our  simulations  show  that  the  elastic  response  of  the  structure  depends 
on  the  ratio  of  the  strut  length  to  its  thickness  ( l/t  in  Figure  4-4(a)),  or  it  s  slenderness 
ratio. 

For  short  struts  with  small  l/t  ratio  (l/t  <  1.67),  the  load  transfer  is  as  follows. 
Figures  4-5(a)  and  (b)  show  the  atomic-level  tensile  and  shear  strain  distribution  for 
the  (45  A,  78  A.  30  A)  nanoporous  structure.  Also,  Figures  4-5(c)  and  4-5(d)  shows 
a  ^-section  slice  of  the  structure  while  undergoing  elastic  deformation.  Both  the  load 
distribution  and  the  deformation  profile  of  the  Z-plane  cross-section,  show  that  the 
vertical  struts  are  mainly  under  pure  tensile  load,  whereas  the  horizontal  struts  arc 
under  pure  shear  load.  The  total  tensile  deformation  in  the  structure  can  thus  be 
written  as, 


we\2  =  (Z  —  t)  sin  <?  +  we22\  (4.1) 

where  t  is  the  thickness  of  the  struts,  w  is  the  height  of  the  vertical  strut  and  l  is 
the  length  of  the  horizontal  strut,  cp  is  the  angle  the  length  of  horizontal  strut  makes 
with  the  A"  axis,  as  shown  in  Figure  4-6a,  and  s22  and  e22  are  the  tensile  strains  in 
the  overall  structure  and  the  vertical  struts  respectively. 

For  small  deformations,  sirup  «  tan  (p  ~  e[2:  whereA12  is  the  shear  strain  in  the 
horizontal  strut.  Thus  Equation  4.1  can  be  written  as: 


ur  _ l—t  l  ,  w  _ l—t  t  i  a  . 

-2  -  £i2  t  %2  -  yy  777  +  yyp 


<JT  i  —  t  T 
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u 22 
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£■22  ' 


(4.2) 


where  r  and  /r12  are  the  shear  stress  and  shear  modulus  for  the  horizontal  strut,  a 
and  E2o  are  the  tensile  stress  and  elastic  modulus  in  the  Y  direction  for  the  vertical 
strut,  and  aT  and  E22  are  the  tensile  stress  and  elastic  modulus  in  the  Y  direction  for 
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Figure  4-5.  Strain  plots  and  deformation  shape  analysis  for  a  nano-honeycomb  structure 
with  thick  struts.  For  the  (45  A.  78  A,  30  A)  structure  at  an  overall  tensile  strain  of  0.105. 
(a)  and  (b)  show  a  measure  of  the  per-atom  tensile  strain  and  shear  strain  through  the 
components  of  the  Jacobian  of  the  local  deformation  matrix,  i.e.,  Jy-y-  and  J\y  respectively. 
Note  that  only  Si  atoms  are  visualized.  The  strain  distributions  clearly  show  the  presence  of 
simple  shear  loading  in  the  horizontal  struts,  and  tensile  loading  in  the  vertical  struts.  For 
the  (45  A,  78  A.  30  A)  structure,  (c)  shows  a  Z  cross-section  of  the  undeformed  structure: 
(d)  shows  the  same  cross-section  at  a  overall  strain  of  0.105.  The  atomic  bonds  clearly  show 
the  simple  shear  loading  in  the  horizontal  struts. 
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Figure  4-6:  Load  distribution  and  stress  transfer  for  the  nano-honeycombs  with  thick 
struts;  (a)  shows  the  various  geometry  parameters  used  to  calculate  elastic  strains  in  the 
honeycomb  structure,  and  the  deformed  shape  for  structures  with  thick  struts,  and  (b) 
shows  the  load  balance  between  the  vertical  and  the  horizontal  struts. 
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the  overall  structure.  Doing  a  force  balance  on  the  members,  as  shown  in  Figure  4- 
6(b).  gives  us,  at  —  2rt  =k  a  —  2 r.  And,  a  force  balance  overall  gives  us,  at  —  2aTl. 
Combining  these  with  Equation  4.2  gives  us  the  elastic  modulus  of  the  nanoporous 
material, 


i  21  n-t  i  i  \  (ui-t)\  i  /2 nj_ 

Coo  t  V  re  2pi2  +  E22)  V  wt  )  ^12  +  \  t  )  E22' 

To  check  the  validity  of  Equation  (4.3),  by  treating  the  geometry  parameters, 
1(1  —  t)  /wt,  and  21 /t  as  variables  we  perform  a  least  square  curve  fit  over  four  different 
nano-honeycomb  geometries  for  the  overall  structure  elastic  modulus  obtained  from 
molecular  dynamics  simulations  (Table  4.1).  These  give  us  a  shear  modulus  of  the 
struts  as  13.3  GPa  and  elastic  modulus  as  109.1  GPa.  These  can  be  compared  to 
bulk  a-quartz  ReaxFF  values  of  shear  modulus  of  14.1  GPa,  Young’s  modulus  (E22) 
of  67.4  GPa  and  elastic  constant,  C22  of  130.3  GPa.  The  calculated  E22  lies  between 
the  Young’s  modulus  (Poisson  effect  in  the  transverse  directions)  and  the  clastic 
constant  C2 2  (no  Poisson  contraction)  showing  that  the  vertical  struts  have  boundary 
conditions  between  these  two  extremes. 


honeycomb 

slenderness 
ratio  (l/t) 

1(1  —  t)/(tw) 

2  l/t 

E(GPa) 

(45  A,  78  A,  30  A) 

1.38 

0.29 

2.75 

20.06 

(67  A,  112  A,  30  A) 

1.33 

0.30 

2.67 

22.5 

(56  A,  78  A,  30  A) 

1.20 

0.15 

2.40 

30.38 

(67  A,  89  A.  30  A) 

1.17 

0.13 

2.33 

31.98 

Table  4.1:  Geometric  parameters  and  clastic  moduli  for  nano-honeycomb  structures  with 
thick  struts  (l ft  <  1.67).  For  the  nomenclature  for  referring  to  the  structure  specifications, 
refer  to  Figure  4-4(a). 


The  pores  change  shape  as  the  structure  is  deformed.  The  change  in  shape  can 
be  followed  by  the  changing  angle  between  a  vertical  strut  and  an  adjacent  horizontal 
strut,  with  strain.  The  elastic  modulus  equation  (Equation  (4.3))  can  also  be  checked 
by  measuring  the  angle  of  deviation  of  the  horizontal  struts  from  their  original  hor¬ 
izontal  position  as  a  function  of  strain.  Equation  4.3  predicts  the  following  relation 
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between  the  sine  of  the  angle  and  applied  strain. 


£22  ~  ( -  +  j  yin  0-  (4.4) 

V  w  {E22/V 12)/ 

This  plot  is  shown  in  Figure  4-7.  for  a  particular  nano-honeycomb  structure,  the 
(45  A,  78  A.  30  A).  Interestingly  the  sine  of  the  angle  variation  remains  linear  with 
applied  strain,  as  predicted  from  Equation  4.4  during  elastic  deformation,  but  the 
angle  becomes  constant  on  the  onset  of  plasticity  in  the  form  of  shearing  in  the 
struts. 


Figure  4-7:  The  variation  in  pore  shape,  measured  through  sine  of  the  angle  between 
the  length  of  the  horizontal  strut  and  the  X  axis  (sec  Figure  4-6).  as  a  function  of  applied 
strain.  The  relation  is  linear  on  the  regime  of  elastic  deformation  (part  of  the  curve  in 
diamond-shaped  data  points  colored  blue),  as  predicted  by  the  theory.  There  is  no  change 
in  the  angle  during  plastic  deformation  (part  of  the  curve  in  square  data  points  colored 
pink).  This  shows  that  the  pore  shapes  remain  about  the  same  during  plastic  deformation. 


For  structures  with  more  slender  struts,  with  larger  l/t  ratio  ( l/t  >  1.67),  the 
response  of  the  horizontal  struts  are  seen  to  be  governed  by  bending  deformation. 
Figure  4-8(a)  and  (b)  show  Z-section  slice  of  the  structure  while  undergoing  elastic 
deformation  for  the  (45  A,  133  A,  30  A)  nanoporous  structure.  The  deformation 
profile  of  the  Z-cross-section  show  that  the  vertical  struts  are  mainly  under  pure 
tensile  load,  whereas  the  horizontal  struts  are  under  beam  bending  load. 

In  this  case,  the  horizontal  struts  can  be  approximated  as  beams  fixed  at  both 
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Figure  4-8:  Deformation  shape  analysis  for  a  nano-honeycomb  structure  with  slender 
struts.  For  the  (45  A,  133  A,  30  A)  structure,  (a)  shows  a  Z  cross-section  of  the  undeformed 
structure;  (b)  shows  the  same  cross-section  at  an  overall  tensile  strain  of  0.133.  The  atomic 
bonds  clearly  show  the  bending  loading  in  the  horizontal  struts. 


ends,  with  a  point  load  in  the  middle.  The  horizontal  displacement  can,  in  this  case 
be  written  as: 


A  _  P(2/):i 

192*» /!  3  (4.5) 

where  P  =  a  t ;  I  =  jp . 

In  this  equation.  A  is  the  vertical  displacement  at  the  end  of  a  horizontal  strut 
(see  Figure  4-9).  P  is  the  load  on  the  beam,  1  is  the  moment  of  inertia  of  the  strut 
(assuming  its  width  in  the  Z  direction  =1),  and  Ell  is  the  elastic  modulus  of  the 
horizontal  strut  in  the  X  direction.  This  simplifies  to: 


A  = 


aP 


2  Ent°-' 


The  total  tensile  deformation  in  the  structure  can  thus  be  written  as, 


(4-0) 


ICffoo  = 


a 


P 


2  EP 


+  W£%o- 


Or.  the  total  strain  in  the  structure, 


(4.7) 
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Figure  4-9:  Load  distribution  and  stress  transfer  for  the  nano-lioncycombs  with  slender 
struts;  (a)  shows  the  various  geometry  parameters  used  to  calculate  elastic  strains  in  the 
honeycomb  structure,  and  the  deformed  shape  for  slender  struts,  and  (b)  shows  the  load 
balance  and  deflection  relation  between  the  vertical  and  the  horizontal  struts. 
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(4.8) 


_T 
’  22 


(7  <7 

2EuWt2  E22 


Thus  the  elastic  modulus  can  be  calculated  as, 


1  /  Z*  \  1  /2/\  1 

El  ~  \wt*J  Eu  +  [t)  E22  (4'9) 

To  check  the  validity  of  Equation  (4.9).  treating  the  geometry  parameters,  lA/wt 3, 
and  2 l/t,  as  variables,  we  perform  a  least  square  curve  fit  over  three  different  nano- 
honcycomb  geometries  with  slender  struts  for  the  overall  structure  elastic  modulus 
obtained  from  molecular  dynamics  simulations  (Table  4.2).  These  give  us  an  Eu 
of  the  struts  as  76.2  GPa  and  E22  as  104.3  GPa.  These  can  be  compared  against 
ReaxFF  values  of  bulk  a- quartz  Young’s  modulus  in  the  X  direction  (Eu)  of  65.8 
GPa  and  Young’s  modulus  in  Y  direction  (£22)  of  67.4  GPa  and  elastic  constant 
C22  of  130.3  GPa.  The  calculated  E2 2  here  too  lies  between  the  Young’s  modulus 
(Poisson  effect  in  the  transverse  directions)  and  the  elastic  constant  C22  (no  Poisson 
contraction)  showing  that  the  vertical  struts  have  boundary  conditions  between  these 
two  extremes.  The  calculated  Eu  lies  close  to  the  Young’s  modulus  in  the  X  direction, 
showing  that  the  horizontal  struts  under  bending  load  arc  almost  completely  under 
free  surface  boundary  conditions. 


honeycomb 

slenderness 
ratio  (l/t) 

Z4/ (wt3) 

2  l/t 

E(GPa) 

(17  A,  78  A.  30  A) 

2.83 

21.48 

5.67 

2.97 

(45  A,  134  A,  30  A) 

2.00 

9.14 

4.00 

6.41 

(33  A.  78  A.  30  A) 

1.67 

3.86 

3.33 

11.89 

Table  4.2:  Geometric  parameters  and  elastic  moduli  for  nano-honeycomb  structures  with 
slender  struts  (l/t  >  1.67).  For  the  nomenclature  for  referring  to  the  structure  specifications, 
refer  to  Figure  4-4(a). 
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4.6.2  Plasticity  and  failure 


All  nano-honevcomb  structures  fail  either  by  brittle  crack  propagation  or  plastic  de¬ 
formation  followed  by  necking  in  the  struts.  Here,  we  find  the  critical  size  parameters 
at  which  we  observe  this  brittle- to-ductile  transition.  We  observe  that  as  the  width 
of  the  constituent  strut  elements  (parameter  t  in  Figure  4-4a)  is  reduced,  the  large 
strain  deformation  behavior  changes  from  crack  initiation  and  propagation  to  plas¬ 
tic  deformation  in  some  of  the  struts.  This  change  is  seen  for  strut  sizes  of  ps  60 
A,  and  the  change  is  gradual,  with  some  intermediate  nano-honeycomb  structures 
showing  presence  of  both  plasticity  and  crack  propagation  in  different  areas.  The 
structures  exhibiting  plastic  deformation  show  crystalline  shearing  of  thin  struts.  We 
also  observe  that  once  plasticity  starts  the  flow  stress  does  not  change  until  one  of 
the  deforming  struts  starts  to  form  a  neck  region.  Once  a  strut  ruptures,  whether 
through  brittle  or  ductile  fracture,  other  struts  still  carry  load  and  thus  the  overall 
structure  shows  stress  softening  and  graceful  failure  (as  seen  in  the  last  part  of  the 
stress-strain  curves  in  Figure  4-2). 

The  plastic  deformation  of  the  struts  is  seen  to  occur  due  to  crystal  slip,  which 
will  happen  at  the  theoretical  shear  strength.  This  is  in  conformity  with  the  idea 
of  a  flaw-tolerant  size  in  brittle  materials,  where  below  a  certain  crystal  size  the 
presence  of  cracks  does  not  affect  the  fracture  stress  and  the  material  always  fails  at 
t  he  theoretical  strength  [33].  Using  this  and  the  stress-transfer  path  in  the  structure, 
we  can  predict  the  yield  strength. 


a 


r 


yield  A h  j  ■ 


(4.10) 


where  <Jyield  is  the  yield  stress  of  the  structure  and  rt/Js  the  theoretical  shear 
strength  in  silica.  Figure  4-10  shows  a  fit  of  the  yield  stress  versus  the  structure 
geometry  (t/l),  which  allows  a  prediction  of  the  theoretical  shear  strength,  which  for 
silica  turns  out  be  ~  5.46  GPa.  To  calculate  a  length  scale  order-of-magnitude  at 
which  the  flaw-tolerant  behavior  steps  in,  we  have  to  find  a  strut  size  below  which 
plastic  deformation  sets  in.  Assuming  a  maximum  stress  concentration  at  a  crack 
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tip  of  a  semi-infinite  crack  in  an  infinite  solid  of  height  t,  and  since  the  plastically 
deforming  strut  is  under  shear  loading,  using  mode  II  loading  equations  on  the  crack, 
we  get  through  application  of  the  Griffith  criterion, 


G=Tr2^  {i-n) 

where  G  is  the  energy  release  rate  at  the  crack  tip,  andys  is  the  fracture  surface 
energy.  So  to  induce  plasticity,  the  system  stress  should  reach  the  theoretical  shear- 
strength,  or, 


t  < 


4p75 


(4.12) 


~  5  -i 


CD 

>  0 


0  0.2  0.4  0.6  0.8 

t/l  (A/A) 


Figure  4-10:  Plot  of  yield  stress  versus  a  geometry  parameter  for  nano-honeycomb  struc¬ 
tures  showing  plastic  flow.  Equation  (4.10)  in  the  main  text  shows  that  the  slope  of  the 
linear  fit,  here  5.46  GPa,  corresponds  to  the  theoretical  shear  strength  of  the  struts. 


Substituting  values  for  the  (0001)  cleaved  surface  energy  in  silica  (0.17eV/A2) 
[155],  and  the  shear  modulus  and  the  theoretical  shear  strength  (from  Equation  4.10). 
we  get  an  estimate  of  ~  80  A.  This  is  close  to  the  strut  size  of  60  A  below  which 
we  see  evidence  of  plasticity  in  our  molecular  dynamics  simulations.  Thus  for  t  strut 
thickness  values  less  than  this  value,  the  nano- honeycomb  structures  will  exhibit 
plasticity. 
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4.7  Conclusions 


The  mechanics  of  diatom-inspired  nanoporous  silica  structures  under  tensile  loading 
has  been  investigated  using  molecular  dynamics  simulation  and  theoretical  analysis. 
We  find  that  the  elastic  modulus  of  the  nanoporous  structures  is  geometry  dependent, 
and  can  be  modified  for  a  given  porosity  by  changing  the  microstructure.  Under  elastic 
deformation,  the  load  distribution  in  the  structure  is  dependent  on  the  slenderness 
ratio  of  the  constituent  struts.  For  thick  struts  (l/t  <  1.67),  the  structure  carries 
stress  through  a  tension-shear-tension  loading  chain  where  the  horizontal  struts  are 
under  pure  shear  loading  and  the  vertical  struts  under  tensile  load.  For  slender 
struts  (l/t  >  1.67),  the  structure  passes  load  through  a  tcnsion-bencling-tcnsion  load¬ 
ing  chain;  the  vertical  struts  are  again  under  tensile  load  while  the  horizontal  struts 
behave  like  beams  with  fixed  ends  under  bending  load. 

The  structure  changes  shape  as  it  is  loaded,  with  the  pore  shapes  deforming.  The 
pore  shape  change  can  be  tracked  by  an  angle  change  between  adjoining  struts,  and 
it;  is  shown  that  there  exists  a  linear  relation  between  the  sine  of  this  angle  0  (Figure 
4-7(a)),  and  the  applied  strain  in  the  regime  of  elastic  deformation. 

On  further  loading,  the  structures  either  undergo  plastic  deformation  followed  by 
necking  or  brittle  fracture  starting  by  crack  propagation  from  corners  of  the  pores. 
There  exists  a  minimum  strut  width  size  below  which  plasticity  can  be  seen  («  60  —  80 
A),  and  this  size  is  justified  based  on  the  flaw-tolerance  concept  for  brittle  materials 
failure.  For  structures  with  strut  widths  larger  than  this,  fracture  always  proceeds  by 
brittle  crack  initiation  and  propagation.  For  strut  widths  below  this  size,  plasticity 
is  seen  in  the  form  of  shear  deformation  in  the  horizontal  struts  under  shear  loading. 
The  yield  strength  of  these  plastically  deforming  structures  has  been  calculated  to 
depend  on  the  structure  geometry. 

An  important  consideration  for  the  mechanical  response  of  nanoporous  silica  struc¬ 
tures  is  the  effect  of  presence  of  water,  in  particular  for  a  more  in-depth  consideration 
of  the  properties  of  diatoms  as  they  live  in  aqueous  environments.  Preliminary  in¬ 
vestigations  by  Garcia  et  al.  [156]  suggest  that  nanoporous  silica  structures  (similar 
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to  those  found  in  this  study)  have  reduced  ductility  and  yield  stress  once  the  surface 
is  hydrated.  Similar  results  have  been  obtained  in  a  previous  study  of  silica  nanorod 
deformation  in  the  presence  of  water  using  semi-empirical  quantum  mechanics  meth¬ 
ods  [157].  The  authors  concluded  that  strained  siloxane  fSi-O-Si)  bonds  arc  attacked 
by  water  which  results  in  lower  stress  and  lower  failure  strain  of  the  silica  nanorod, 
compared  to  a  dry  silica  nanorod. 

Our  results  establish  a  size-dependent  brittle-to-ductile  transition  in  nanoporous 
silica  structures.  Similar  behavior  has  been  observed  experimentally  in  silica  nanowires, 
where  the  size  range  for  the  transition  is  20  nm.  The  aspect  ratio  and  shape  of  pores 
can  be  modified  to  change  the  yield  strength  of  these  structures.  The  high  values 
of  stress  for  yielding  or  fracture  lead  to  large  enhanced  ductility  in  these  materials 
over  bulk  silica.  The  structures  that  show  plastic  yielding  also  show  large  toughness 
improvement  over  bulk  silica.  These  results  reveal  that  nano-scaling  with  control  of 
porous  geometry  can  lead  to  application  of  silica  in  carrying  loads  in  small  devices. 
The  increased  toughness,  elastic  ductility  and  plastic  ductility  arising  from  nano¬ 
scaling  may  also  be  fundamental  in  understanding  the  use  of  similar  structures  by 
nature  in  creating  porous  exoskeletons  in  diatoms. 

Apart  from  the  improved  ductility,  one  property  that  is  affected  to  a  large  extent 
is  the  material  stiffness.  Nanostructuring,  as  outlined  in  this  chapter,  can  reduce  the 
stiffness  of  the  silica,  structures  by  up  to  90%  of  its  original  value.  Is  it  possible  to 
design  structures  which  recover  the  stiffness  value  of  bulk  silica  while  retaining  the 
toughness  improvement  of  the  nanoporous  silica?  In  the  next  chapter,  we  design  hi¬ 
erarchical  silica  nanocompositcs  with  2  levels  of  hierarchy.  Wc  proceed  by  developing 
an  mesoscale  method  for  studying  the  mechanics  of  silica  structures  at  the  micron 
length  scale.  We  then  use  the  method  in  the  analysis  of  the  mechanics  of  hierarchical 
silica  structures. 
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Chapter  5 


Mesoscale  Model  of  Deformation 
and  Failure  of  Hierarchical  Silica 


N  anocomposites 


This  chapter  is  focused  on  the  development  of  a  mesoscale  model  for  studying  the 
mechanics  of  hierarchical  silica  nanocomposite  structures.  Biomineralized  silica-based 
materials,  such  as  sea  sponge  exoskeletons,  and  diatoms,  showcase  the  use  of  inferior 
constituent,  materials  such  as  brittle  ceramics  like  silica  and  soft  protein  to  create 
materials  that  have  surprisingly  high  toughness  and  resistance  to  crack  propagation, 
and  retain  stiffness  values  close  to  the  ceramic  constituent  (sec  Chapter  1).  A  common 
design  paradigm  seen  in  those  materials  is  the  existence  of  multiple  levels  of  structural 
hierarchy  [22,  29,  23.  28,  31].  Various  studies  have  pointed  at  how  the  existence  of 
these  structural  design  levels  and  their  functional  adaptation  helps  the  material  in 
retaining  the  best  possible  combination  of  properties  of  the  constituent  materials 
[2,  158]. 


5.1  Review  of  structural  bio-silica  materials 

A  fascinating  example  of  biomineralized  structures  are  those  of  diatoms  [140,  159, 
160],  a  microscopic  algae  that  feature  cell  walls,  or  frustules,  mainly  composed  of 
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amorphous  silica.  A  porous  hierarchical  structure  dominates  the  landscape  of  the  cell 
wall  and  encompasses  intricate  patterns  that  are  highly  varied  and  ordered,  as  shown 
in  Figure  1-2,  and  reach  all  the  way  down  to  the  nanoscale.  Please  see  Chapter  1, 
section  1.2  and  1.3  for  more  details  on  the  diatom  structure  and  mechanical  properties. 
Another  silica-based  design  is  found  in  the  deep  sea  sponge  ( Euplectella  sp.)  skeleton 
(Figure  5-1).  The  structure  shows  structural  design  at  several  length  scales,  ranging 
from  a  composite  consisting  of  consolidated  nanometer  sized  silica  spheres  embedded 
in  an  organic  matrix  at  the  sub-micron  scale  (Figure  5- 1(d)),  through  cylindrical 
lamellar  silica-organic  material  composite  at  a  micron  length  scale  (called  spicules- 
Figure  5- 1(c))  to  a  square-lattice  cage-like  structure  consisting  of  cylindrical  rods 
(Figure  5-l(a.-b))  at  the  macroscale. 


Figure  5-1:  Structural  hierarchies  in  a  silica-based  skeletal  structures  in  a  sea  sponge,  (a) 
shows  the  external  cage  structure  of  the  silica-based  skeletal  system  of  Euplectella  sp.,  scale 
bar  5  mm  (b-d)  show  some  of  the  underlying  hierarchical  structures  with  (b)  showing  fiber- 
composite  structure  in  a  constituent  beam  consisting  of  many  spicules,  scale  bar  20  microns, 
(c)  single  spicule  showing  laminated  silica-protein  structure,  scale  bar  5  microns  and  (d) 
biosilica  constituent  of  the  silica  layers  revealing  its  consolidated  nanoparticlc  nature,  scale 
bar  500  nm.  Figure  reproduced  from  [28]. 


We  summarize  some  key  mechanical  properties  of  diatoms  and  sea  sponge  ex¬ 
oskeletons  here.  Several  studies  reported  in  the  recent  literature  have  revealed  the 
mechanical  properties  ol  diatom  shells.  Hamm  et  al.  [39]  used  a  glass  needle  to  load 
and  break  diatom  frustulcs  in  order  to  probe  their  mechanical  response  at  failure, 
and  found  high  strength  and  reversible  clastic  strains  (e.g.  2.5%  reversible  strain  in  a 
frustule  section).  Other  researchers  [40]  have  used  AFM  nanoindentation  to  study  the 
nanoscale  material  properties  of  the  porous  frustule  layers  of  diatoms,  identifying  pore 
sizes  on  the  order  of  several  tens  of  nanometers  at  the  smallest  levels  in  the  hierarchy, 
with  ultra-thin  silica  walls  on  the  order  of  several  nanometers.  They  observed  that 
the  variation  of  mechanical  properties  between  the  hierarchical  frustule  layers  could 
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be  influenced  by  the  pore  size,  pore  distance,  porosity,  and  under  different,  biominer¬ 
alization  processes.  Experimental  tests  to  identify  the  mechanical  properties  of  sea 
sponge  spicules  have  also  been  reported.  The  spicules  consisting  of  97-98%  volume 
fraction  of  silica  and  2-3%  of  organic  material,  arc  seen  to  be  tougher  in  both  tension 
and  bending  than  silica  glass,  by  a  factor  of  6-7  times  [161].  The  thin  organic  layers 
present  between  silica  layers  in  these  spicules  show  delamination,  crack-bridging,  and 
provide  elastic,  viscoelastic  and  viscoplastic  means  for  energy  absorption  under  load¬ 
ing,  leading  to  the  high  toughness  [28,  161].  The  fracture  strength  was  also  improved 
3-5  times  [162,  163].  with  large  improvement  in  ductility  over  silica  glass  [163.  30]. 

We  hypothesize  that  the  nanoporons  geometry  and  hierarchical  arrangement  of 
the  frustulcs  in  diatoms,  and  the  hierarchical  structure  of  the  spicules  and  their  ar¬ 
rangement  in  sea  sponge  skeletons  are  crucial  to  providing  enhanced  toughness  at  high 
strength  and  stiffness  even  though  the  constituting  material  itself  (that  is.  silica)  is 
inherently  brittle  and  mechanically  inferior  for  structural  applications.  Previous  work 
by  the  authors  has  focused  on  nanoporous  structures  of  silica  and  shown  that  this 
geometry  leads  to  high  ductility  and  toughness  [151,  164].  although  at  the  cost  of 
structural  stiffness.  We  thus  explore  here  the  creation  of  a  multi-level  hierarchical 
composite  structures  made  of  nanoporous  silica  and  much  stiffer  bulk  silica,  and  in¬ 
vestigate  their  mechanical  properties  through  an  in  silico  approach. 

The  focus  of  the  study  reported  in  this  article  is  first  the  development  of  an 
atomisticallv-informed  mesoscale  particle-spring  model,  which  is  then  used  for  mod¬ 
eling  these  composite  structures  at  the  micro-scale  [165].  The  method  is  tested  on 
the  mechanics  of  randomly-dispersed  fibcr-roinforccd  composite  structures  made  from 
bulk  silica  and  nanoporous  silica  phases,  to  enable  fracture  studies  at  micrometer 
length-scales.  An  important  question  our  model  is  expected  to  answer  is  whether  or 
not  it  is  possible  to  create  a  highly  functional  material  with  great  strength  and  tough¬ 
ness  out  of  a  single  material  constituent.,  silica,  by  solely  engineering  its  structural 
design  at  multiple  levels  without  adding  any  additional  material  component. 


Ill 


5.2  Materials  and  Methods 


A  multiscale  bottom-up  computational  methodology  is  used  here  to  study  the  effect 
of  hierarchical  design  on  material  properties  and  the  mechanics  of  deformation.  At 
the  nanoscale  we  use  molecular  dynamics  simulations  to  study  the  mechanics  of  the 
nanoporous  silica  structures.  Molecular  dynamics  with  the  first  principles  based  reac¬ 
tive  ReaxFF  atomistic  force  field  is  a  powerful  tool  to  capture  fundamental  nanoscale 
phenomena  and  the  mechanisms  behind  them.  At  the  micron  length  scale,  we  de¬ 
velop  a  mcsoscalc  spring-lattice  network  model.  The  model  is  derived  from  as  well  as 
validated  against  the  atomistic  results.  Spring-lattice  network  models  at  the  micron 
length  scales  arc  able  to  capture  elasticity,  plasticity  and  fracture  phenomena  at  these 
length  scales.  We  describe  the  details  of  the  methods  in  the  following  sections. 

At  the  nanoscale,  the  first-principles  derived  ReaxFF  force  field  [62]  is  used  to 
characterize  the  mechanical  behavior  of  nanoporous  silica  structures,  as  discussed  in 
Chapter  4.  Fully-atomistic  simulations  have  been  carried  out  for  the  mechanics  of 
nanoporous  silica  structures  [151,  164]  (sec  Chapter  4).  These  studies  show  the  effect 
of  pore  size,  distribution  and  porosity  on  elastic  modulus,  plasticity,  ductility  and 
toughness  of  the  structures.  Figure  5-2(a)  shows  one  of  these  characteristic  na no¬ 
honeycomb  silica  structures.  The  availability  of  these  studies  allows  us  to  extract 
constitutive  laws  of  nanoporous  silica  behavior  that  can  be  used  to  build  mesoscale 
models  of  silica  structures  with  hierarchies. 

5.2.1  Mesoscale  method  development  and  validation 

Figure  5-2(b)  shows  the  model  setup  consisting  of  a  network  of  material  particles 
connected  in  a  lattice  arrangement  through  springs.  Such  two-dimensional  spring- 
lattice  networks  have  been  used  previously  to  model  deformation  and  fracture  in 
brittle  and  quasi-brittle  materials  [166,  167,  168,  169],  and  are  particularly  suitable 
for  studying  fracture  phenomena  in  heterogeneous  materials  [170,  171,  172],  The 
two-dimensional  nature  of  the  model  resembles  plane  strain  loading  conditions.  The 
constitutive  stress-strain  law  under  tensile  load  is  obtained  for  a  particular  nano- 
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Figure  5-2:  Atomistic  and  mesoscale  modeling  approaches  are  combined  here  to  describe 
the  material  from  the  nano-  to  the  micro-scale.  Parameters  for  the  mesoscale  model  are 
derived  from  constitutive  behavior  at  the  nanoscale  obtained  using  atomistic  simulations. 
Panel  (a)  shows  the  geometry  of  the  nano-honeycomb  used  as  building  blocks  for  the  com¬ 
posite  structures,  panel  (b)  shows  a  section  of  the  triangular  mesh  mesoscale  particle-spring 
model  setup,  panel  (c)  show  stress-strain  curves  obtained  from  atomistic  simulations  of  a 
nano-honeycomb  structure,  and  for  bulk  silica  with  a  crack  of  the  same  size  as  the  pores 
in  the  nano- honeycomb.  The  legend  defines  the  nano- honeycomb  structure,  which  is  shown 
as  (t,  pi,  pw )  parameters  for  the  structure  (numerical  values  given  in  A).  The  bulk  silica 
structure  shows  purely  brittle  fracture,  the  nano-honeycomb  structure  show  ductile  frac- 
ture.  Panel  (d)  shows  the  behavior  of  the  mesoscale  triangular  mesh  lattice  fitted  to  this 
constitutive  behavior  (the  agreement  with  the  full  atomistic  result  depicted  in  panel  (c)  is 
evident  and  provides  direct  validation  of  the  mesoscale  model). 


honeycomb  silica  and  bulk  silica  using  atomistic  simulations,  shown  in  Figure  5-2(c). 
In  this  figure,  the  bulk  silica  structure  has  a  pre-crack  with  dimensions  of  the  pore 
size  in  the  nano- honeycomb  structure,  to  compare  structures  with  similar  defect  sizes. 

As  seen  in  Figure  5-2(c).  the  bulk  structure  is  stiff  and  brittle,  while  the  nanoporous 
structure  is  soft  and  ductile.  The  force-extension  law  for  the  mesoscale  inter-particle 
potential  is  hyperelastic  and  is  fit  to  the  constitutive  law  behavior  of  nano-silica  and 
bulk  silica  under  tensile  load  (Figure  5-2(d)).  The  hyperelastic  spring  potential  mod¬ 
els  the  atomistic  results  for  the  nano-honeycomb  as  elastic-perfectly  plastic  behavior, 
and  the  bulk  silica  as  elastic-brittle  behavior.  For  the  nano-honeycomb,  the  flow  stress 
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is  obtained  from  atomistic  simulations,  and  is  calculated  as  the  mean  stress  during 
plastic  deformation.  Since  the  aim  of  this  study  is  to  obtain  mechanics  of  different 
composite  structures,  the  properties  of  the  local  springs  in  the  lattice  in  a  certain 
material  model  are  changed  according  to  whether  they  lie  geometrically  inside  the 
matrix  or  reinforcing  phase. 

The  brittle  bulk  silica  phase  inter-particle  potential  is  modeled  as  follows: 

f  kAALA  if  ALa  <  ALAc. 

Fa(ALa)={  (5.1) 

^  0  if  ALa  >  AL4X. 

where  FA  is  the  force  on  a  spring  between  two  material  particles  of  the  brittle 
phase,  kA  is  the  force  constant  for  the  spring,  A LA  is  the  extension  in  the  spring, 
and  A LA  c  is  the  critical  cutoff  distance  for  the  spring.  The  cutoff  distance  ALa_c, 
beyond  which  the  mesoscale  spring  carries  zero  load,  is  fixed  based  on  the  failure 
strain,  and  the  force  constant  kA  is  fit  to  the  elastic  modulus  of  the  material  obtained 
from  atomistic  simulations. 

The  ductile  nano-porous  phase  inter-particle  potential  is  modeled  as: 


Fb  (AT)  - 


kB  AL  B  if  A Ljg  *£  A Lb.i  * 

kBAL  b.  i  +  kc  (A  Lb  —  ALbA)  if  ALb.i  <  A  Lb  <  A  Lb.c  « 

0  if  ALb  >  ALb,c  ■ 


(5.2) 

where  FB  is  the  force  on  a  spring  between  two  material  particles  of  the  ductile 
phase.  ALb  is  the  extension  in  the  spring,  and  A LBA  is  the  extension  for  the  spring 
for  the  onset  of  the  plastic  regime,  and  ALb,c  is  the  critical  cutoff  extension  when  the 
spring  breaks  and  stops  carrying  load,  kB  is  the  force  constant  for  the  elastic  response, 
and  kc  =  kBj  100  is  a  small  force  constant  to  model  plastic  deformation  at  constant 
flow  stress.  A LBA  is  tit  to  the  yield  strain,  and  A Lb.c  is  fit  to  the  failure  strain,  and 
kB  is  fit  to  the  elastic  modulus  for  the  nano-honeycomb  structure  obtained  from  the 
atomistic  simulations  (Figure  5-2(c)). 

The  mesoscale  model  uses  a  triangular  lattice  regular  mesh  for  the  arrangement 


114 


of  the  springs.  This  result s  in  isotropic  elasticity  behavior  which  can  be  assumed 
if  the  underlying  structure  is  polycrystalline  with  regular  arrangement  of  grains. 
Anisotropic  behavior  could  also  be  introduced  in  the  mesoscale  model  in  principle, 
but  for  the  proof-of-conccpt  study  targeted  in  this  study  (in  the  spirit  of  simple  com¬ 
putational  experiments),  we  assume  isotropic  behavior  in  the  underlying  atomistic 
structure.  The  implications  of  this  assumption  are  that  crystal-orientation  dependent 
anisotropic  elastic  and  fracture  behavior  phenomena  will  not  be  captured  through  this 
model.  The  interfaces  between  the  two  phases  in  the  composite,  if  the  crystal  struc¬ 
ture  is  not  continuous  across  the  interface,  would  also  possibly  contribute  to  plasticity 
in  the  material  through  slip  and  friction,  however  the  interfaces  have  been  just  as¬ 
signed  the  low  strength  and  stiffness  values  of  the  nanoporous  phase  in  the  model. 
This  doesn’t  allow  the  model  to  capture  any  toughness  enhancement  by  interfacial 
plasticity  n rechanisms . 

Bulk  two-dimensional  mesoscale  models  are  constructed  of  bulk  silica  and  nanoporous 
silica  and  subject  to  tensile  testing.  Comparison  of  the  elastic  moduli  and  fracture 
toughness  between  the  atomistic  simulations  and  mesoscale  simulations  are  used  to 
fix  spring  constants  and  inter-particle  distance  in  the  mesoscale  model. 

The  clastic  modulus  and  fracture  toughness  for  the  atomistic  model  are  calcu¬ 
lated  for  a  bulk  silica  sample  with  a  ccntcr-orack.  The  fracture  toughness  using  the 
ReaxFF  force  held  for  silica  is  calculated  to  be  0.79  MPa-yhn.  This  is  rather  close  to 
experimental  values  in  the  literature  for  fused  quartz,  0.6-0.75  MPa-v/m  [173,  174], 
The  elastic  modulus  for  the  spring-lattice  model  is  fit  to  102.3  GPa.  and  the  mode  I 
fracture  toughness  to  be  0.79  MPayhn. 

To  match  the  atomistic  simulation  values  (see  brief  review  above),  an  inter-particle 
distance  of  78  nm  and  spring  constants  of  3,932  N/m  and  134.4  N/m  for  the  brittle 
and  ductile  phases,  respectively,  arc  chosen  for  a  through  thickness  of  100  nm.  This 
ensures  for  a  separation  of  scales  between  the  scales  described  by  the  atomistic  model 
and  the  characteristic  length-scale  associated  with  the  mesoscale  model  (the  typical 
scales  of  the  atomistic-level  models  is  up  to  ten  nanometers).  All  mesoscale  models 
are  implemented  within  the  LAMM  PS  software  package  [175]. 
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It  is  noted  here  that  the  model  parameters  can  easily  be  adapted  to  describe 
other  nanostructures,  and  can  even  be  extended  to  describe  multiple  nanostructures 
in  the  study  of  hierarchical  systems.  The  specific  model  considered  here  represents 
one  specific  case  study  explored  here. 


5.2.2  Fracture  property  characterization 

For  materials  that  show  failure  by  growth  of  a  single  dominant  crack,  we  characterize 
them  by  calculating  their  fracture  toughness.  The  toughness  is  calculated  by  calcu¬ 
lating  the  energy  release  rate  by  invoking  the  J-integral  [176,  177]  in  its  domain  form 
[178]  around  the  crack  tip,  given  by: 


J  = 


dut  \ 

dx2) 


dq 

dX~ 


dS , 


(5.8) 


where  W  is  the  strain  energy  density,  Pr)  is  the  first  Piola-Kirchoff  stress  tensor. 
Ui  represents  the  displacement  field.  X  arc  the  material  coordinates,  Sq  represents 
the  undeformed  area  of  the  domain,  5  is  the  Kroneeker  delta,  and 


( 1  = 


0  :  /•  /'|. 

1  :  r  =  r2, 

(r  —  )  /  (r2  —  7*i )  :  7~i  <  r  <  r2. 


(5.4) 


where  the  parameters  rx  and  ro  are  shown  in  Figure  5-9 (a).  The  discrete  form  of 
this  equation,  for  small  displacements,  is  given  by  (see,  e.g.  [179,  180]): 


2  r 

J=EE 

aeSn  i,j  L 


was2j 


pa  dUj  \  dq(Xa 

13  dx. 


dX3 


So, 


(5.5) 


where  Sft  is  the  initial  undeformed  area  occupied  by  the  material  particle  a ,  Xa 
is  the  initial  position  of  material  particle  a,  Wa  is  the  local  strain  energy  density  at 
any  material  particle  a  which  is  calculated  as  follows, 


whore  <pa  is  the  potential  strain  energy  of  the  material  particle  a,  obtained  from 
simulation  by  splitting  the  spring  potential  energy  between  material  particles  sharing 
the  spring  bond,  and  the  stress  at  the  atom  a.  ( P°j )  is  calculated  from  the  virial 
theorem  [181]: 

p"*wEr"'s®f,S’  <5'7) 

where  rQ'3  is  a  vector  joining  material  particles  a-  and  ,5,  ia3  is  the  force  applied  on 
material  particle  a  by  material  particle  8,  and  =  Sat  is  the  volume  occupied  by 
material  particle  a,  where  Sa  is  the  area  occupied  by  the  material  particle  a  in  the 
deformed  configuration.  This  formulation  of  the  J-integral  is  used  to  avoid  involving 
high  stress  values  at  the  crack  tip  region  in  the  calculation,  and  the  convergence  of 
the  J-integral  is  checked  by  measuring  its  value  against  different  integration  domain 
regions.  The  strain  and  duL/dX\  are  obtained  by  a  local  least  square  fit  to  the 
neighbor  displacement  field  at  each  material  particle  location. 

5.2.3  R-curve  calculation 

Stable  crack  advance  for  every  load  configuration  is  noted  by  finding  the  crack  tip 
location.  A  particular  spring  bond  is  regarded  as  broken  when  its  deformation  exceeds 
the  cutoff  for  the  interaction.  Crack  surfaces  are  visualized  by  finding  all  spring  bonds 
which  have  snapped  for  a  given  load.  The  J-integral  is  used  to  find  the  energy  release 
rate  for  a  given  amount  of  stable  crack  advance.  Plot  of  the  J-integral  from  the 
start  of  crack  initiation  through  crack  propagation  provides  the  R-curve  [182]  for  the 
material,  i.e.  how  its  fracture  toughness  changes  as  a  function  of  stable  crack  advance. 

5.3  Results  and  discussion 

We  consider  models  of  randomly-distributed  fiber-reinforced  composites  of  bulk  silica 
with  small  volume  fractions  of  nanoporous  silica.  Figure  5-3  shows  the  different 
geometries  considered  here.  In  all  cases,  the  fibers  are  circular  in  cross-section  and 
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randomly  distributed  through  the  matrix,  oriented  parallel  out-of-plane  (in  the  Z 
direction  in  Figure  5-3)  and  loaded  in  the  cross-section  (X-Y  plane).  Fiber  diameters 
are  1.3  pm  in  a  model  size  of  27  pm  by  23.3  pm  and  100  nm  out-of-plane.  Two 
kinds  of  composites  arc  designed-  a),  hard-fibers  of  bulk  silica  (of  volume  fraction 
76%)  embedded  in  a  soft  nanoporous  matrix,  and  b),  soft  fibers  of  nanoporous  silica 
dispersed  in  a  hard  silica  matrix  (of  volume  fraction  86%).  The  two  designs  arc  shown 
in  Figure  5-3.  We  select  five  different  random  structures  for  both  morphologies  as 
representative  bulk-silica  rich  composites  to  obtain  statistically  relevant  data,  and  will 
refer  to  them  henceforth  as  the  brittle-fiber  and  brittle-matrix  structures,  respectively. 


Z<$ - >X 


Figure  5-3:  Geometry  of  randomly  distributed  fiber-composite  structures  at  the 

mcsoscalc.  Constituents  arc,  bulk  silica  (in  grey/light,  with  a  high  volume  fraction)  and 
nano-honeycomb  structures  (in  blue/dark,  small  volume  fraction).  Design  conditions  that 
enhance  toughness  of  bulk  silica  by  distribution  of  small  amount  of  nano-honeycomb  silica 
are  being  investigated  here.  The  structure  on  the  left  shows  bulk  silica  fibers  as  reinforce¬ 
ment.  whereas  the  right  structure  shows  nano-honeycomb  silica  fibers  as  the  reinforcing 
phase.  The  structure  sizes  arc  27  pm  by  23.3  pm  and  100  nm  out-of-plane.  Reinforcing 
phase  is  always  in  fiber  form  of  diameter  1.3  pm,  aligned  parallel  out-of-plane  but  randomly 
distributed  in-planc.  Structures  of  both  types  arc  studied  for  crack  propagation  response 
under  mode  I  loading  with  plane  strain  conditions.  Initial  pre-crack  sizes  range  from  3.9  to 
5.4  pm.). 
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We  now  present  a  systematic  study  of  the  effect  of  composition  and  material 
distribution  on  the  elasticity,  plasticity  and  fracture  of  these  composites.  Figure  5-4 
shows  the  elastic  modulus  measurements  for  both  the  brittle-fiber  and  brittle-matrix 
materials.  We  observe  that  for  a  random  dispersion  of  the  second  phase  in  the  first, 
the  elastic  modulus  lies  between  the  Rcuss  and  Voigt  bounds.  However,  the  brittle- 
matrix  composite  is  seen  to  have  a  modulus  closer  to  the  upper  bound  (Voigt)  than 
the  brittle-fiber  composite.  Thus  brittle-matrix  morphologies  are  more  suitable  for 
obtaining  a  higher  stiffness  for  random  fiber  distributions. 


Figure  5-4:  Elastic  modulus  for  the  brittle-fiber  and  brittle-matrix  composites  plotted  as 
a  function  of  volume  fraction  of  the  brittle  phase.  The  limits  at  0  and  1  volume  fraction 
correspond  to  nano-honeycomb  silica  (soft)  and  bulk  silica  (stiff)  respectively.  The  com¬ 
posite  moduli  are  seen  to  lie  within  the  Reuss  and  Voigt  bounds  of  load  sharing  between 
the  phases,  the  brittle-matrix  morphology  providing  stiffness  closer  to  the  upper  limit  (of 
Voigt  modulus). 


Next  we  create  sharp  edge  cracks  in  all  materials  and  load  them  under  quasi¬ 
static  mode  I  loading.  Loading  is  carried  out  by  stepped  edge  displacement  boundary 
conditions  and  relaxing  t  he  global  positions  of  all  material  particles  using  a  conjugate 
gradient  energy  minimization  scheme  [183],  Initial  crack  size  is  3. 9-5.4  pm.  Crack 
initiation  is  identified  bv  the  advance  of  the  crack  front  at  a  particular  loading  strain 
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value.  Plastic  deformation  is  visualized  by  calculating  local  strains  in  springs  and 
comparing  to  yield  strain  values.  We  notice  the  absence  of  any  plasticity  prior  to  crack 
propagation  initiation  in  all  the  composite  designs.  Figure  5-5  shows  representative 
overall  stress-strain  response  for  composite  structures  with  brittle-fiber  and  brittle- 
matrix  morphologies,  with  and  without  the  presence  of  the  edge  crack. 


strain 


Figure  5-5:  Stress-strain  curves  for  (a)  composite  structures  with  bulk-silica  as  the  rein¬ 
forcement,  with  and  without  presence  of  a  pre-crack.  The  near-identical  response  shows  the 
flaw-tolcrance  behavior  for  these  structures  to  pre-cracks  of  the  given  size.  The  structures 
show  multiple  cracking  throughout  the  material,  and  this  is  reflected  in  the  stress-strain 
curve  as  a  gradual  loss  of  stiffness  of  the  material  as  the  number  and  size  of  the  multiple 
cracks  grow,  (b)  Stress-strain  curves  for  composite  structures  with  nano-honeycomb  sil¬ 
ica  as  the  reinforcement,  with  and  without  presence  of  a  pre-crack.  The  varying  fracture 
strengths  clearly  show  an  effect  of  the  crack  size.  All  structures  fail  by  the  growth  of  a 
single  dominant  crack. 


The  brittle-fiber  structures  show  multiple  micro-cracking  sites  (throughout  the 
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Figure  5-6:  Crack  pathways  for  composite  structures  with  nano-honeycomb  structure  as 
the  matrix  and  brittle  silica  as  the  reinforcing  fiber  phase.  The  volume  fraction  of  silica 
phase  is  76%.  (a)  and  (b)  show  fracture  progress  starting  from  a  material  with  no  pre-crack; 
(c)  and  (d)  show  fracture  progress  in  the  same  material  with  a  pre-crack  present.  In  these 
cases  we  observe  that  the  pre-crack  propagates  for  a  small  distance  but  does  not  propagate 
through  the  sample,  and  other  smaller  cracks  are  initiated  throughout  the  sample.  These 
multiple  small  cracks  determine  the  stress-strain  response  of  the  structure.  The  structure 
is  thus  flaw-tolerant  to  pre-cracks  of  these  sizes,  and  the  fracture  stress  and  behavior  are 
almost  independent  of  the  size  of  the  pre-crack.  This  is  reflected  also  in  the  stress-strain 
curve  for  the  stress-strain  response  of  the  undefected,  and  cracked  structures  (shown  in 
Figure  5-5a).  The  fracture  toughness  cannot  be  measured  for  these  structures  and  crack 
sizes. 


sample)  under  tensile  load.  The  presence  of  pre-cracks  is  not  seen  to  affect  this  phe¬ 
nomenon,  and  the  material  fails  by  the  growth  and  coalescence  of  several  micro-cracks 
(Figure  5-6).  This  phenomenon  is  also  captured  by  the  absence  of  any  effect  on  the 
stress-strain  curve  to  the  presence  of  pre-cracks  of  a  certain  size  and  below  (Figure 
5-6(a)),  and  can  be  classified  as  a  defect- tolerant  state.  Fracture  mechanics  formula¬ 
tions  cannot  be  used  for  such  material  microstructures  and  crack  sizes,  and  damage 
mechanics  which  deals  with  evolution  of  damage  with  applied  load  for  example  in 
the  form  of  diffuse  micro-cracking  would  have  to  be  used  to  characterize  the  failure 
response.  1  he  diffuse  cracking  also  leads  to  an  effect  on  the  elastic  modulus  (slope 
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of  the  stress-strain  curve)  of  the  material,  as  seen  in  Figure  5-6(a),  the  modulus  is 
lowered  at.  higher  strain  values  when  the  microcracks  start  to  increase  in  size  and 
number. 
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Figure  5-7:  Crack  pathways  for  composite  structures  with  brittle  silica  as  the  matrix  and 
nano-honeycomb  structure  as  the  reinforcing  fiber  phase.  The  volume  fraction  of  bulk  silica 
Phase  is  86%.  (a)  and  (b)  show  fracture  progress  starting  from  a  no  pre-crack  material; 
(c)  and  (d)  show  fracture  progress  from  the  same  material  with  a  pre-crack.  In  both 
cases  wc  observe  that  fracture  occurs  through  the  propagation  of  a  dominant  crack,  (a) 
and  (c)  show  the  un-cracked  and  pre-cracked  specimens  at  the  same  load,  the  un-cracked 
specimen  is  intact,  whereas  the  pre-crack  has  started  propagating  in  the  other  specimen. 
Since  the  fracture  strength  of  a  structure  with  a  dominant  propagating  crack  is  pre-crack- 
size  dependent  (according  to  fracture  mechanics),  the  stress-strain  curve  for  the  stress-strain 
response  of  the  undcfcctcd,  and  cracked  structures  arc  markedly  different  (shown  in  Figure 
5-5b).  Fracture  toughness  can  be  measured  for  these  structures,  as  the  energy  required  for 
the  growth  of  the  pre-crack  per  unit  crack  advance 


The  brittle-matrix  structures,  on  the  other  hand,  show  failure  by  growth  of  a 
single  dominant  crack  (Figure  5-7).  This  phenomenon  is  captured  in  the  stress-strain 
curve  by  the  decrease  in  fracture  stress  with  increase  in  pre-crack  size  (Figure  5- 
5(b)),  as  expected  from  fracture  mechanics.  Snapshots  of  the  pre-crack  growth  are 
shown  in  Figure  5-8  for  all  the  brittle-matrix  composite  structures  considered.  The 
crack  path  shows  several  phenomena,  such  as  linking  fibers  lying  in  its  path,  and 
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following  a  tortuous  pathway.  In  some  cases,  we  see  crack  branching  at  a  fiber-matrix 
interface  and  in  one  case  (Figure  5-8(b))  the  crack  is  not  continuous  but  is  bridged  by 
a  ligament  of  the  matrix.  These  highlight  all  the  different  mechanisms  at  this  length 
scale  that  can  affect  the  fracture  behavior  of  the  material. 


0.48%  strain  0.48%  strain  0.48%  strain 


0.48%  strain  0.48%  strain 

Figure  5-8:  Different  composite  structures  with  brittle  bulk-silica  as  the  matrix  (volume 
fraction  86%)  and  ductile  nano-honeycomb  structures  as  the  reinforcing  fiber  phase  showing 
fracture  toughness  improvement  mechanisms.  The  fibers  have  circular  cross-section  and  are 
randomly  distributed  and  five  different  random  structures  arc  shown  here.  A  single  pre-crack 
is  introduced  and  then  subjected  to  mode  I  loading.  Propagation  of  the  single  dominant 
crack  is  seen  on  loading,  and  the  propagation  path  is  marked  in  white.  All  structures  show 
that  the  crack  path  is  not  straight,  but  connects  reinforcing  fibers  lying  close  to  the  original 
crack  plane.  Crack  deflection  and  bridging  by  reinforcing  fibers  behind  the  crack  tip  are 
the  mechanisms  seen  to  increase  toughness  here. 

The  fracture  toughness  is  calculated  for  the  brittle-matrix  composite  structures 
using  the  J-intcgral  formulation  (Figure  5-9(a)).  Crack  initiation  toughness  is  mea¬ 
sured  by  a  J-integral  calculation  of  the  strain  energy  release  rate  just  prior  to  crack 
propagation.  R-enrve  measurements  (crack  advance  resistance  as  a  function  of  stable 
crack  growth)  are  then  undertaken  by  obtaining  stable  crack  propagation  distances 
for  different  load  strains  and  repeating  the  j-integral  calculations.  We  thus  obtain 
the  fracture  toughness  as  a  function  of  crack  advance  distance.  We  notice  all  struc- 
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Figure  5-9:  Calculation  ol  the  J-integral  and  R-curves.  (a)  shows  the  J-integral  calculation 
for  a  stationary  crack  by  the  use  of  a  domain-integral  around  the  crack.  The  J-integral 
provides  the  value  of  the  energy  release  rate  per  unit  advance  of  the  crack  into  the  crack 
tip.  oi  the  lcsistancc  to  crack  propagation.  The  red/dark  region  shown  is  the  domain  of 
integration  and  the  convergence  of  the  J-integral  is  tested  by  taking  different  and  regions 
for  the  same  crack  and  specimen  configuration. (b)  shows  fracture  toughness  measures  as 
a  function  of  crack  advance  (R-curve  behavior)  for  all  the  structures  in  Figure  5-8.  The 
toughness  of  bulk  silica  is  also  shown  as  a  dotted  curve. 


turcs  show  an  initiation  toughness  close  to  bulk  silica  and  significant  improvement 
in  toughness  («  4.4  times)  as  the  crack  proceeds  (Figure  5-9(b)).  The  improve¬ 
ment  in  toughness  arises  from  erack-deficction  and  crack-bridging  by  the  reinforcing 
phase,  and  also  ligament-bridging  by  the  matrix.  Simultaneously,  the  stiffness  is  not 
compromised  due  to  the  large  volume  fraction  of  the  bulk  silica  phase  (Figure  5-4). 
Steady-state  fracture  propagation  toughness  values  can  be  obtained  from  the  R-curve 
for  the  composite  for  large  crack  growths  in  the  range  of  several  microns. 
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5.4  Summary  and  Conclusions 


We  established  here  a  hierarchical  mcsoscalc  modeling  scheme  that  represents  the 
mechanics  of  silica-nanoporous  silica  composites  at  the  micron  length  scale,  bridging 
scales  from  a  few  nanometer  (accessible  to  full- atomistic  modeling)  to  scales  of  tens  of 
micrometers.  The  particle-spring  mesoscale  model  has  parameters  fitted  to  atomistic 
simulations  using  t  he  ReaxFF  force  held.  The  ReaxFF  force  held  itself  has  parameters 
fitted  to  first-principles  quantum  simulations  of  energetics  of  silicon-oxygen  clusters. 
We  establish  the  validity  of  the  mesoscale  model  by  capturing  the  stress-strain  relation 
of  the  bulk  silica  and  nanoporous  silica,  and  the  fracture  toughness  and  fracture 
mechanisms  (e.g.  crack  bridging,  spreading,  etc.)  of  silica  in  different  hierarchical 
designs. 

Enabled  by  the  mesoscale  model,  we  studied  the  fracture  mechanics  of  randomly- 
dispersed  fiber  composite  structures  using  this  coarse-grained  model.  ReaxFF  simula¬ 
tions  have  shown  us  previously  that  the  design  of  nanoporous  silica  leads  to  enhanced 
ductility  and  toughness  at  the  cost  of  compromising  stiffness.  We  demonstrated  that 
through  using  composites  with  a  large  volume  fraction  of  the  bulk  silica,  it  is  possible 
to  obtain  structures  with  stiffness  almost  as  large  as  the  silica  phase,  but  signifi¬ 
cant  toughness  improvement  due  to  the  use  of  nanoporous  materials  as  the  other 
phase.  To  this  end.  we  completed  a  case  study  of  fracture  of  two  kinds  of  compos¬ 
ites.  one  characterized  by  hard  fibers  of  bulk  silica  embedded  in  a  soft  nanoporous 
matrix  (brittle-fiber),  and  another  one  that  consists  of  soft  fibers  of  nanoporous  silica 
dispersed  in  a  hard  silica  matrix  (brittle-matrix). 

The  brittle-fiber  structures  showed  multiple  micro-cracking  initiation  under  tensile 
load,  and  the  material  fails  by  their  growth  and  coalescence,  independent  of  presence 
of  pre-cracks  in  the  material  (Figure  5-6).  Fracture  mechanics  formulations  cannot  be 
used  for  such  material  microstructures  and  crack  sizes,  and  damage  mechanics  which 
deals  with  evolution  of  damage  with  applied  load  e.g.  in  the  form  of  diffuse  micro¬ 
cracking,  have  to  be  used  to  characterize  the  failure  response  in  this  flaw-tolerant 
regime  of  cracking. 
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The  brittle-matrix  structures,  on  the  other  hand,  showed  failure  by  growth  of  a 
single  dominant  crack  (Figure  5-7).  The  fracture  toughness  was  calculated  for  several 
randomly  dispersed  fiber-composite  structures  using  the  J-integral  formulation.  The 
crack  toughness  is  tracked  as  a  function  of  crack  advance,  using  an  R-curve  (Figure 
5-9 (b)),  and  we  observed  that  all  structures  show  an  initiation  toughness  close  to 
bulk  silica  and  significant  improvement  in  toughness  as  the  crack  proceeds,  similar  to 
findings  reported  earlier  for  bone  [184],  The  improvement  in  toughness  arises  from 
crack-deflection  and  crack-bridging  by  the  reinforcing  phase,  and  ligament  bridging 
by  unbroken  matrix.  At  the  same  time,  the  stiffness  of  the  composite  structure  lies 
close  to  its  Voigt  upper  limit  (Figure  5-4).  To  the  best  of  our  knowledge,  this  type 
of  R-curve  analysis  is  reported  in  this  article  for  the  first  time  for  hierarchical  silica 
structures,  and  implemented  in  a  multiscale  simulation  scheme  of  these  structures.  We 
note  that  our  model  is  simple  and  has  several  limitations,  such  as  the  lack  of  describing 
anisotropies,  and  its  two-dimensional  nature.  Nevertheless,  it  provides  a  useful  tool 
for  computational  experiments  geared  to  elucidate  fundamental  design  principles  of 
biological  materials  at  multiple  length-scales,  and  the  formulation  presented  here  can 
be  easily  adapted  to  other  cases. 

The  improvement  in  toughness  and  retaining  of  stiffness  of  bulk  silica  by  designing 
small  regions  of  nanoporous  geometry  in  the  bulk  phase  point  to  a  design  method¬ 
ology  for  obtaining  stiff  and  tough  materials  out  of  an  inherently  brittle  material 
(silica)  to  begin  with  (Figure  5-10).  The  design  philosophy  can  be  summed  up  as 
the  use  of  a  single  material  silica  which  is  traditionally  considered  an  undesirable 
(mechanically  inferior)  structural  material  due  to  its  brittleness,  arranged  in  a  hier¬ 
archical  pattern  with  substructures  that  go  down  to  the  nanoscale  dimensions.  A 
future  refinement  in  this  model  would  be  modeling  mechanical  response  of  interfaces 
between  the  phases  in  the  composite  using  atomistic  simulations  and  incorporating 
interfacial  parameters  in  the  mesoscale  model.  Further  work,  shown  in  Chapter  6, 
aims  at  studying  effect  of  additional  levels  of  hierarchies  on  these  properties.  The 
present  design  scheme  suggested  enables  us  to  obtain  highly  functional  materials  that 
feature  enhanced  toughness  and  stiffness.  The  utilization  of  the  design  paradigm 
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Figure  5-10:  The  design  of  the  hierarchical  silica  composite  improves  the  toughness  of 
bulk  silica  significantly  (~  4.4  times)  while  compromising  on  the  stiffness  only  slightly  (~  70 
%  of  bulk).  This  points  towards  the  use  of  hierarchies  along  with  a  single  design  material 
to  improve  undesirable  mechanical  properties  significantly  (here  low  toughness)  while  not 
compromising  on  the  desirable  ones  (here  high  stiffness). 
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outlined  here  leads  to  a  substantial  increase  the  design  space  for  brittle  materials, 
without  the  need  to  add  additional  materials,  and  solely  by  geometrical  design. 

Likely  driven  by  evolutionary  pressures,  materials  in  living  systems  already  use 
this  design  philosophy  in  a  variety  of  protective  armor  and  internal  load-carrying 
structures  in  many  different  species  through  the  use  of  hierarchies  (as  it  was  shown 
lor  soft  protein  materials  in  earlier  studies  [185]),  and  our  simulations  point  to  the 
strength  of  using  this  philosophy  for  designing  engineering  structures  made  out  of 
hard  materials  as  well.  The  exploitation  of  this  method  for  de  novo  material  design 
could  be  substantial,  as  it  provides  a  pathway  to  utilize  nanoscale  material  elements 
with  superior  properties  for  macroscale  highly  functional  materials — with  superior 
properties,  despite  the  utilization  of  simple  and  abundant  building  blocks  such  as 
silica.  In  the  next  C  hapter,  we  extend  the  use  of  the  mesoscale  model  developed 
heie,  to  the  study  of  fracture  properties  and  toughness  of  hierarchical  materials  with 
multiple  levels  of  hierarchy.  Our  attempt  is  to  observe  numerically,  how  the  addition 
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of  more  levels  of  hierarchy  improves  toughness  properties. 
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Chapter  6 


The  role  of  multiple  structural 
hierarchy  levels  in  defect  tolerance 
and  toughness 


This  chapter  extends  the  use  of  the  mesoscale  model  developed  in  Chapter  5  to  the 
direct  simulation  of  the  fracture  properties  of  structures  with  several  levels  of  hierar¬ 
chy.  This  is  motivated  by  the  observation  that  many  biological  structural  materials 
show  not  just  2  or  3^  hut  several  levels  of  structural  hierarchies  from  the  nanoscale 
up  to  the  macroscale.  Cortical  bone,  for  example  has  7  distinct  structural  hierarchy 
levels  [1,  186],  using  two  basic  constituents  of  collagen  protein  and  hydroxyapatite 
mineral.  Sea  sponge  exoskeletons,  made  of  silica  and  protein,  show  6  levels  of  struc¬ 
tural  hierarchy  [28].  A  question  that  arises  is  how  does  the  increase  in  number  of 
hierarchy  levels  impact  mechanical  properties?  Does  it  give  rise  to  properties  that 
are  not  apparent  from  1  or  2  levels  of  hierarchy?  Here,  we  will  specifically  focus  on 
the  affect  of  adding  more  hierarchical  levels  on  crack  advance  and  toughness. 
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6.1  Background  on  structures  with  multiple  lev¬ 
els  of  structural  hierarchy  and  study  of  their 
mechanics 

In  biological  structural  materials  possessing  multiple  levels  of  structural  hierarchies, 
many  phenomena  and  properties  are  linked  across  structural  scales  (see  Chapter  1. 
section  1.3).  Experimentally,  fracture  behavior  of  these  materials,  such  as  bone  and 
nacre  have  been  shown  to  be  heavily  linked  across  scales  [35,  30],  with  mechanisms  at 
several  different  hierarchy  levels  participating  in  the  overall  behavior.  This  has  also 
led  to  some  contention  as  to  which  mechanisms  are  the  important  ones  [184,  33]. 

The  field  of  computational  and  theoretical  studies  of  mechanics  of  materials  with 
a  large  number  of  hierarchical  levels  is  in  its  infancy.  A  major  problem  in  the  field 
is  the  issue  of  studying  deformation  mechanisms  and  their  interactions  across  a  huge 
range  of  length  scales,  which  are  often  beyond  the  computational  capabilities  of  a  sin¬ 
gle  computational  technique  such  as  atomistic  simulation  or  finite- clement  methods. 
Another  challenge  is  that  a  representative  volume  element  cannot  be  found  at  any 
length-scale  without  accounting  for  its  entire  sub-structure  down  to  the  nanoscale, 
because  of  the  heterogeneity  of  the  design  at  any  length-scale.  This  heterogeneity 
can  potentially  lead  to  very  large  models  in  terms  of  computational  cost,  beyond  the 
reach  of  today's  computational  resources.  Multi-scale  methods  (as  described  in  Chap¬ 
ter  2)  have  to  be  applied  in  these  cases.  Nonetheless,  concurrent  multi-scale  methods 
cannot  be  applied  easily  due  to  the  non-applicability  of  localized  deformation;  and 
hierarchical  multi-scale  methods,  tend  to  develop  into  too  coarse  a  description  at  very 
large  length  scales,  due  to  the  loss  of  detailed  information  about  mechanisms  with 
every  additional  level  of  coarse-graining. 

One  recent  approach  has  been  the  study  of  self-similar  hierarchical  assemblies 
with  some  assumptions  about  similarity  in  failure  mechanisms  at  each  hierarchy  level 
[22.  30,  187,  188].  These  continuum  models  predict  strength,  stiffness  and  toughness 
scaling  with  number  of  hierarchy  levels  for  a  number  of  simple,  self-similar  assemblies 
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e.g.  for  self-similar  hierarchical  assembly  of  bone  nanostructure  [36]  or  gecko  adhe¬ 
sion  [189.  190].  A  key  drawback  of  these  models  is  the  use  of  a,  self-similar  geometry 
assumption,  and  also  the  assumption  of  correspondence  of  energy  dissipation  mecha¬ 
nisms  and  failure  mechanisms  across  hierarchies.  Self-similarity  in  geometric  design 
is  not  seen  in  many  biological  materials,  such  as  bone,  nacre,  sea  sponge  exoskeleton, 
diatoms  etc.  [186,  158,  28,  27],  where  the  structural  arrangement  at  each  hierarchy 
level  is  very  different.  Another  approach  has  been  of  taking  into  account  explicit 
probabilistic  modes  of  failure  pathways  down  several  levels  of  hierarchy,  assuming 
certain  possible  unit  failure  events  (e.g.  failure  of  unit  hydrogen  bonds  in  hierarchical 
alpha- helix  structures,  and  their  combinatorial  arrangements)  [191,  192],  Some  qual¬ 
itative  models  have  also  been  derived  to  identify  design  similarities  in  the  hierarchies 
of  various  biological  materials  that  may  offer  the  key  to  improvements  of  properties 
[13.  31,  193], 

In  this  work,  our  aim  is  to  study  the  effects  of  hierarchical  arrangements  on  tough¬ 
ness.  Toughness  in  materials  is  often  measured  through  fracture  toughness,  which  fo¬ 
cuses  on  energy  dissipated  in  the  structure  per  unit  length  advance  of  a  major  crack 
that  eventually  leads  to  complete  failure  of  the  structure.  However,  as  demonstrated 
in  the  last  chapter,  complex  microstruetures  can  have  very  different  failure  mecha¬ 
nisms  than  the  propagation  of  a  single  large  crack,  e.g.  through  diffuse  damage  in 
the  form  of  several  cracks  interacting  throughout  the  structure.  We  aim  to  study 
how  toughness  is  enhanced  through  hierarchies,  specifically  given  that,  through  the 
fracture  mechanics  viewpoint,  the  crack  will  have  to  interact  with  all  different  lev¬ 
els  of  hierarchy.  Specifically,  our  goal  is  to  observe  whether  each  additional  level  of 
structural  hierarchy  effect  has  a  substantial  effect  on  toughness. 

We  will  also  study  a  related  property,  the  defect-tolerance  length  scale.  Defect 
tolerance  measures  the  sensitivity  of  fracture  strength  to  the  presence  and  size  of  a 
crack.  A  higher  defect-tolerance  implies  a  lower  sensitivity  to  crack  size,  i.e.,  a  large 
change  in  the  size  of  a  crack  present  in  the  material/structure  is  required,  for  a  small 
change  in  the  fracture  strength.  One  of  the  aims  of  robust  design  is  the  increase  of  the 
defect-tolerance  length  scale,  so  that  the  material/structure  strength  does  not  vary 
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for  a  range  of  flaw  size  distribution.  This  can  be  related  to  the  parameter  insensitivity 
definition  of  robustness  in  biological  systems  theory  [25]  (see  Chapter  1,  section  1.2). 
Through  the  use  of  flaw-tolerant  design  at  the  nanoscale  in  brittle  materials  [33] .  and 
suppressing  dislocation-mediated  plasticity  in  ductile  materials  [129],  it  is  possible  to 
eliminate  any  dependence  of  fracture  strength  on  crack  size  for  very  small  nanomctcr- 
sized  crystals.  Here,  however,  we  also  want  to  study  the  effect  of  hierarchies  on  this 
sensitivity  of  fracture  strength  on  crack  size. 

With  this  aim.  we  use  the  bulk-silica / nanoporous  silica  model  developed  in  Chap¬ 
ters  4  and  5.  We  focus  on  the  use  of  silica  as  a  ‘poor'  structural  material,  and  attempt 
to  engineer  it  to  possess  better  properties,  through  nanostructuring  (Chapter  4)  and 
use  of  hierarchies  (Chapter  5  and  here).  In  this  chapter,  the  atoinistiftally- informed 
mesoscale  method  developed  in  Chapter  5,  is  thus  extended  to  several  levels  of  hier¬ 
archy. 

6.2  Mesoscale  simulation  results 

In  many  biological  materials  possessing  several  hierarchy  levels,  the  design  substruc¬ 
ture  at  each  level  is  periodic  and  built  from  an  repeating  template.  Examples  include 
cortical  bone,  nacre  at  all  their  hierarchy  levels  (see  Chapter  1,  section  1.2);  diatoms 
and  sea-sponge  skeletons  at  all  hierarchy  levels  except  the  nauoscale  biosilica,  which 
is  a  randomly-dispersed  composite  (see  Figures  1-2  and  5-1). 

6.2.1  Two-hierarchy  level  structures  with  periodic  geometry 

We  extend  the  mesoscale  designs  of  Chapter  5.  to  regularly  distributed  (periodic) 
composite  structures.  The  protein  constituent  found  in  biological  structures,  which 
is  soft  and  tough  is  replaced  by  nanoporous  silica  in  our  structures;  the  mineral 
constituent,  which  is  hard  and  brittle  is  replaced  by  bulk  silica.  For  a  particle- 
reinforced  composite,  there  are  two  distinct  design  schemes  possible,  one  in  which 
the  soft  material  is  the  matrix,  the  other  where  it  is  the  reinforcement  particles.  We 
thus  choose  two  representative  systems,  one  in  which  the  soft/tough  material  is  the 
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continuous  phase,  with  the  hard/brittle  material  dispersed  in  platelet  form  within 
the  matrix,  mimicking  the  periodic  arrangement  observed  in  several  hierarchy  levels 
in  bone  and  nacre  [1].  The  second  representative  system  is  where  the  hard/brittle 
constituent  is  the  continuous  phase,  with  regions  of  soft/tough  material  embedded 
within  the  hard  matrix.  These  structures  arc  found  in  biocrystals,  where  protein 
material  is  encapsulated  in  a  hard  crystal,  such  as  biological  calcite  single  crystals 
[194,  195].  The  use  of  a  periodic  distribution  of  one  material  in  another  not  only 
allows  us  to  capture  actual  design  morphologies  observed  in  biological  structures,  but 
also  removes  the  stochastic  element  of  the  dependence  of  mechanical  properties  on 
random  morphologies. 

These  regularly  distributed  composite  structures  arc  shown  in  Figure  6-1.  The 
volume  fraction  of  the  stiff  silica  phase  is  kept  at  a  high  value  and  constant  at  80%. 
For  the  bone-like  arrangement,  the  platelets  of  bulk  silica  have  a  rectangular  shaped 
cross-section  and  8.4  pm  by  2.4  pm  in  size  in  the  X-Y  plane,  and  100  nm  in  the 
out-of-plane  direction.  They  are  arranged  in  a  staggered  fashion,  with  an  overlap 
equal  to  half  their  length  across  subsequent  layers.  The  overall  structure  size  is  27 
pm  by  70  pm  and  100  mil  out-of-plane.  For  the  biocalcite-like  arrangement,  the  soft 
nanoporous  silica  is  embedded  as  rectangular  inclusions  with  a  cross-section  size  of  8.7 
pm  by  0.7  pm  in  size  in  the  X-Y  plane,  and  100  nm  in  the  out-of-plane  Z-direction. 

Next  we  create  sharp  edge  cracks  in  all  structures  and  load  them  under  quasi¬ 
static  mode  1  loading.  Loading  is  carried  out  by  stepped  edge  displacement  boundary 
conditions  and  relaxing  the  global  positions  of  all  material  particles  using  a  conjugate 
gradient  energy  minimization  scheme  [183].  Initial  crack  sizes  range  from  5-20  pm. 
Crack  initiation  is  identified  by  the  advance  of  the  crack  front  at  a  particular  applied 
strain.  Figure  6-2  shows  representative  overall  stress-strain  response  for  composite 
structures  with  bone-like  and  biocalcite-like  morphologies,  with  and  without  the  pres¬ 
ence  of  the  edge  crack.  The  data  shows  that  both  structures  have  a  drop  in  fracture 
stress  in  the  presence  of  a  crack.  The  bone-like  structure  has  a  higher  fracture  stress 
but  shows  a  larger  drop  in  strength  when  a  crack  is  introduced,  whereas  the  biocalcite- 
like  structure  has  smaller  fracture  strength,  but  also  a  lower  sensitivity  to  the  presence 
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Figure  6-1:  Periodic  repeating  design  morphologies  are  chosen  for  the  composite  structure 
at  the  2lld  level  of  hierarchy,  (a)  shows  the  two  representative  systems  chosen,  one  in  which 
the  soft/tough  material  is  the  continuous  phase,  with  the  hard/brittle  material  dispersed 
in  platelet  form  within  the  matrix,  mimicking  several  hierarchy  levels  in  bone  and  nacre: 
the  second  representative  system,  in  the  lower  part  of  (a)  panel,  is  where  the  hard/brittle 
constituent  is  the  continuous  phase,  with  regions  of  soft/tough  material  embedded  within 
the  hard  matiix.  These  structures  arc  found  in  biocrystals,  where  protein  material  is  en¬ 
capsulated  in  a  hard  crystal,  such  as  biological  calcite  single  crystals  [194,  195],  These 
geometric  arrangements  arc  reproduced  using  bulk  silica/nanoporous  silica  materials  using 
the  two-dimensional  mesoscale  model  (see  Chapter  5),  as  shown  in  (b).  which  we  will  refer 
to  as  bone-like  and  biocalcite-like  in  the  subsequent  text. 
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of  a  crack.  Thus  the  biocalcite-like  structure  shows  higher  defect-tolerance  at  the  cost 
of  lower  fracture  strength. 


Figure  6-2:  Stress-strain  curves  for  (a)  bone-like  composite  structures,  with  and  without 
piesence  ot  a  pie-ciack.  Fracture  strength  changes  drastically  on  the  introduction  of  a  crack, 
(b)  Stress-strain  curves  for  a  biocalcite-like  composite  structure,  with  and  without  presence 
of  a  pre-crack.  The  sensitivity  to  fracture  strength  is  much  smaller  than  for  the  bone-like 
composite,  though  the  magnitude  of  the  fracture  strength  is  lower. 


Taking  a  closer  look  at  the  fracture  mechanics  and  crack  propagation  in  both 
structures,  Figure  6-3  gives  us  a  clue  as  to  why  the  biocalcite-like  structure  possesses 
higher  defect  tolerance.  For  small  crack  sizes,  failure  in  these  materials  propagates 
through  the  nucleation  of  several  micro-cracks  at  nano-porous  silica/bulk-silica  inter¬ 
faces  located  from  the  original  crack  tip.  This  makes  the  failure  morphology  appear 
quite  similar  in  the  presence  of  smaller  cracks  and  the  absence  of  any,  and  provides 
comparable  fracture  strength  values.  On  the  other  hand,  for  the  bone-like  structures, 
the  pre-crack  always  propagates  and  is  toughened  by  platelets  bridging  the  wake  of 
the  crack  as  it  propagates. 
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(a) 


(b) 


Figure  6-3:  Crack  pathways  (marked  in  red)  for  hone-like  and  biocalcitc-lic  hierarchical 
structures  in  the  presence  of  a  pre-crack,  (a)  shows  that  for  the  bone-like  structures,  the  pre¬ 
crack  propagates  through  the  sample,  but  the  structure  is  toughened  by  platelets  bridging 
the  wake  of  the  crack  as  it  propagates,  (b)  shows  that  for  small  crack  sizes,  failure  in 
the  biocalcite-like  structure  propagates  through  the  nucleation  of  several  micro-cracks  at 
nano-porous  silica/bulk-silica  interfaces  located  far  from  the  original  crack  tip.  The  fracture 
strength  is  reached  when  several  of  these  micro-cracks  link  up  along  with  the  pre-crack  to 
create  a  failure  pathway  through  the  sample.  These  results  show  that  the  different  crack 
propagation  pathways  in  the  two  structures  lead  to  different  defect-tolerance  response. 
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6.2.2  Extension  to  3-  and  4-hierarchy  level  structures 

The  2-hierarchy  structures  in  the  last  section  are  now  extended  to  an  additional  level 
of  hierarchy.  Both  self-similar  geometries,  in  the  vein  of  earlier  studies  [36]*  and 
dissimilar  geometries  are  considered  (sec  Figure  6-4).  Self-similar  (fractal)  geometries 
are  not  found  in  biological  materials,  but  have  been  used  in  previous  literature  to 
build  simple  models  of  hierarchical  structures;  dissimilar  geometries  are  chosen  as 
more  representative  of  how  the  actual  hierarchy  levels  in  bone,  diatoms  etc.  are  very 
different  from  each  other  in  geometry.  The  self-similar  geometry  uses  a  replica  of  the 
bone-like  arrangement  at  two  different  scales  (Figure  6-4(a)).  The  platelets  size  at 
the  2nd  level  of  hierarchy  is  5.95  pm  by  1.16  pm  in  cross-section,  while  at  the  third 
level  is  12.7  pm  by  5.4  pm.  The  overall  sample  size  is  54  pm  by  70  pm  in  the  X-Y 
plane,  and  100  nm  in  the  out-of-plane  Z-direction.  The  dissimilar  geometry  is  using 
the  bio-calcite  template  for  the  2nd  level,  and  bone-like  arrangement  for  the  3rd  level 
of  hierarchy  (Figure  6-4(b)).  Here,  the  platelets  size  at  the  2nd  level  of  hierarchy  is 
2.02  pm  by  0.47  pm  in  cross-section,  while  at  the  third  level  is  12.6  pm  by  5.5  pm. 
Initial  crack  sizes  in  these  models  range  from  f«5  35  pm. 

Figure  6-5  shows  the  stress-strain  plot  going  from  a  2-hierarchy  bone- like  structure 
to  a  3-hierarchv  self-similar  assembly.  The  volume  fraction  of  the  two  constituents, 
bulk  silica  and  nanoporous  silica,  are  kept  constant  at  80%  and  20%  respectively. 
The  move  from  2-hierarchy  to  3-hierarchy  system  shows,  (a)  a  decrease  in  fracture 
strength,  and  (b)  an  increase  in  defect-tolerance.  Figure  6-6  shows  the  stress-strain 
plot  going  from  a  2-hierarchy  biocalcite-like  structure  to  a  3-hierarchy  dissimilar  as¬ 
sembly.  The  volume  fraction  of  the  two  constituents,  bulk  silica  and  nanoporous  silica, 
are  again,  kept  constant  at  80%  and  20%  respectively.  In  this  case,  the  move  from 
2-hierarchy  to  3-hierarchy  system  shows  a  small  change  in  fracture  strength  decrease 
in  fracture  strength,  and  again  an  increase  in  defect-tolerance. 

We  thus  observe  a  pattern  of  increase  in  defect-tolerance  size  scale  with  the  in¬ 
crease  in  number  of  hierarchy  levels.  In  figure  6-7,  we  take  a  closer  look  at  the 
source  of  this  defect-tolerance.  Here,  we  have  taken  several  samples  of  the  dissimi- 
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Figure  6-4:  Geometry  of  3-hierarchy  composite  structures,  both  self-similar  and  dissimilar 
in  design  between  the  2nd  and  3rd  levels  of  hierarchy.  The  models  are  used  to  study  how  the 
assumption  of  self-similarity  across  hierarchies  affects  stress-strain  response  and  toughness 
values,  (a)  shows  the  self-similar  geometry  using  a  replica  of  the  bone-like  arrangement  at 
two  different  scales.  The  platelets  size  at  the  2nd  level  of  hierarchy  is  5.95  pm  by  1.16  pm  in 
cross-section,  while  at  the  3rd  level  is  12.7  pm  by  5.4  pm.  (b)  shows  the  dissimilar  3-hierarchy 
geometry  using  the  bio-calcite  template  for  the  2nd  level,  and  bone-like  arrangement  for  the 
3rd  level  of  hierarchy.  The  platelets  size  at  the  2lld  level  of  hierarchy  is  2.02  pm  by  0.47  pm 
in  cross-section,  while  at  the  3rd  level  is  12.6  pm  by  5.5  pm.  The  overall  sample  size  is  54 
pm  by  70  pm  by  100  nm  in  the  out-of-planc  direction  (both  cases). 
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Figure  6-5:  Stress-strain  curves  for  (a)  2-hierarchy  bone-like  composite  structures,  with 
and  without  presence  of  a  pre-crack;  (b)  3-hierarchy  self-similar  structure  made  of  bone-like 
composite  structure  at  both  the  2nd  and  3rd  levels,  with  and  without  presence  of  a  pre¬ 
crack.  The  sizes  of  the  pre-cracks  arc  given  in  the  figure  legends.  The  sensitivity  of  fracture 
strength  vs.  crack  size  is  much  smaller  for  the  3-hierarchy  material. 
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Figure  6-6:  Stress-strain  curves  for  (a)  2-hierarchy  biocalcite-like  composite  structures, 
with  and  without  presence  of  a  pre-crack;  (b)  3-hierarchy  dissimilar  structure  made  of 
biocalcite-like  composite  structure  at  the  2nd  level  and  bone-like  at  the  3rd  level,  with  and 
without  presence  of  a  pre-crack.  The  sizes  of  the  pre-cracks  arc  given  in  the  figure  legends. 
The  sensitivity  of  fracture  strength  vs.  crack  size  is  smaller  for  the  3-hierarchy  material. 
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lax  3-hierarchy  structure  with  different  crack  sizes,  and  found  the  fracture  strengths 
(Figure  6-7(a)).  We  find  that  with  a  3-fold  increase  in  crack  size  from  6  pm  to  18 
pm,  there  is  only  a  24%  drop  in  fracture  strength.  A  closer  look  at  the  mechanics 
of  the  stress-strain  curves  is  taken  in  Figure  6-7(b)  and  (c).  Here,  we  observe  that 
the  latter  part  of  the  rising  stress  region  before  fracture  consists  of  the  opening  up  of 
micro-cracks  throughout  the  sample.  Once  these  microcracks  start  moving,  and  link 
up  to  the  pre-crack,  there  is  unstable  crack-propagation  leading  to  a  peak  in  stress 
and  the  fracture  strength.  This  effect  can  also  be  measured  through  the  total  new 
surface  area  created  during  the  diffuse  micro-cracking  regime,  and  while  the  main 
crack  doesn’t  move  (Figure  6-7(e)). 

To  observe  whether  this  pattern  of  increase  in  defect-tolerance  proceeds  with  the 
number  of  structural  hierarchies,  we  create  a  model  with  4  levels  of  hierarchy  (Figure 
6-8).  Here  the  2nd  structural  level  is  biocalcite-like,  while  the  3rd  and  4th  levels 
are  bone-like  in  design.  The  overall  volume  fraction  of  the  bulk-silica  constituent  is 
still  maintained  at  80%.  The  overall  sample  size  is  108  pm  by  140  pm  in  the  X-Y 
plane,  and  100  inn  in  the  ^-direction.  The  system  consists  of  ~3  million  particles 
and  ~9  million  bonds.  The  same  mesoscale  potentials  used  for  the  2-hierarchy  level 
structures  is  still  used,  and  no  further  coarse-graining  in  the  potential  is  carried  out. 
Figure  6-8 (b)  shows  the  stress- strain  curves  for  the  4r hierarchy  structure  with  various 
crack  sizes  from  pm  to  ^64  pm.  Almost  no  change  in  fracture  strength  is  seen 
over  this  very  large  change  in  the  size  of  a  crack  present  in  the  structure.  The 
defect-tolerance,  thus  has  improved  substantially  over  that  of  the  2-hierarcliy  and 
3-hierarchy  structures.  Figure  6-8 (c)  shows  the  R-curve  behavior  over  1,  2,  3.  and  4 
levels  of  hierarchy  structures.  The  R-curve  measures  changes  in  fracture  toughness  as 
a  crack  proceeds  through  the  change  in  energy  released  per  unit  length  stable  crack 
advance. 

We  observe  a  distinct  effect  of  the  addition  of  hierarchies  on  the  R-curve  behavior 
of  the  material.  Both  the  absolute  value  of  Gic  versus  crack  advance,  and  the  slope 
of  the  R-curve,  dGic/dAa  increase  with  the  number  of  hierarchies.  The  slope  of  the 
R-curve  has  a  close  relation  to  the  concept  of  flaw-tolerance.  As  seen  in  Figure  6-9. 
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Figure  6-7:  A  closci  look  at  the  source  of  defect  tolerance  in  3-hicrarchy  structures,  (a) 
shows  the  stress-strain  behavior  and  fracture  strengths  for  several  samples  of  the  dissimilar 
3-hierarchy  structure  with  different  crack  sizes.  We  find  that  with  a  300%  increase  in  crack 
size  from  6  pm  to  18  pm,  there  is  only  a  24%  drop  in  fracture  strength,  (b)  and  (c)  show 
that  the  latter  part  of  the  rising  stress  region  before  fracture  consists  of  the  opening  up  of 
micro-cracks  throughout  the  sample.  These  micro-cracks  arc  shown  in  red  in  part  (b).  with 
the  numbers  indicating  overall  strain  values.  Once  these  microcracks  start  moving,  and 
link  up  to  the  pre-crack,  there  is  unstable  crack-propagation  leading  to  a  drop  in  stress  and 
thus,  the  fracture  strength.  The  lower  part  of  (c)  graph  also  shows  that  this  effect  can  also 
be  measured  through  the  total  new  surface  area  created  during  the  diffuse  micro-cracking 
regime,  at  which  time  the  main  crack  remains  stationary. 
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Figure  6-8:  4-hierarchy  structure  morphology,  stress-strain  plots  and  R-curve  behavior, 
(a)  shows  the  2-hierarchy,  3-hierarchy  and  4-hierarchy  structures  for  comparison,  with  the 
4-hierarchy  structure  having  a  2llfl  hierarchical  level  that  is  biocalcite-like.  while  the  3rd  and 
4  levels  arc  bone-like.  In  the  4-hicrarchy  structure,  the  color  scheme  is:  bulk  silica— red. 
nanoporous  silica  green  and  blue,  to  show  the  4  levels  clearly.  The  overall  volume  fraction 
of  the  bulk-silica  constituent  is  still  maintained  at  80%.  The  overall  sample  size  is  108  pm 
by  140  pm  in  the  X-Y  plane,  and  100  nm  in  the  ^-direction,  (b)  shows  the  stress-strain 
curves  for  the  4-hierarchy  structure  with  various  crack  sizes  from  pm  to  ss64  pm.  Almost 
no  change  in  fracture  strength  is  seen  over  this  very  large  change  in  crack  size.  The  defect- 
tolerance  has  thus  increased  substantially  over  2-hierarchy  and  3-hierarchy  structures,  (c) 
shows  the  R-curve  behavior  over  1,  2,  3,  and  4  levels  of  hierarchy  structures.  The  R-curve 
measures  changes  in  fracture  toughness  as  a  crack  propagates  through  the  change  in  energy- 
released  per  unit  length  of  stable  crack  advance. 


143 


ill  a  material  with  a  rising  R-curve,  unstable  crack  advance  only  occurs  when  the 
following  conditions  are  satisfied. 
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where.  Gapp  is  the  applied  energy  release  rate,  GIc  is  the  fracture  toughness  from  the 
R-curve,  and  a  is  the  crack  length.  This,  along  with  the  load/energy-release  relation 
for  an  edge  crack,  implies  that  the  load  at  which  a  certain  crack  size  will  propagate 
unstably  causing  fracture,  can  be  calculated  by  marking  off  the  crack  size  on  the 
negative- X  axis  of  an  R-curve,  and  constructing  the  tangent  to  the  R-curve  passing 
through  this  point.  The  slope  of  this  curve  is  proportional  to  the  load/  fracture 
stress  at  which  this  crack  propagates  unstably.  This  also  implies  that  for  a  rising 
R-curve,  the  higher  the  rate  of  rise  with  crack  advance,  the  lesser  the  sensitivity  of 
fracture  stress  to  crack  size,  and  thus  the  higher  the  defect-tolerance  length  scale. 
Defect- tolerance  length  scales  are  thus  closely  linked  to  not  only  absolute  values  of 
fracture-crack  initiation  values  [36] ,  but  also  to  the  rising  part  of  the  R-curve  through 
the  slope  of  the  R-curve. 

This  relation  is  clarified  through  the  results  shown  in  Figures  6-10  and  6-11.  Frac¬ 
ture  stress  is  measured  for  2-,  3-  ,  and  4-hierarchv  structures  for  different  crack  sizes. 
The  loss  of  strength  as  a  percentage  of  the  strength  of  the  no-cracked  samples  has 
been  plotted  in  Figure  6-10  for  all  the  levels  of  hierarchy.  We  observe  that  the  sensi¬ 
tivity  of  fracture  strength  to  crack  size  goes  down  with  increasing  hierarchy  level.  If 
a  10%  loss  in  strength  is  fixed  as  a  defect-tolerance  length  scale,  Figure  6-11  shows 
how  this  length  scale  increases  non-linearlv  with  the  number  of  hierarchies. 


6.3  Discussion  and  Conclusions 

In  this  study,  we  study  the  fracture  mechanics  of  hierarchical  structures  with  2 
4  levels  of  hierarchy,  using  an  atomistically-infornied  mesoscale  method.  We  span 
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Figure  6-9:  R-curve  behavior  in  a  material  with  a  rising  R-curve  resistance,  and  link  to 
unstable  crack  propagation.  All  hierarchical  structures  shown  in  Figure  6-8(a)  show  rising 
R-curve  behavior,  where  the  toughness  grows  with  crack  propagation  from  its  initiation 
value.  The  load,  at  which  a  certain  crack  size  will  propagate  unstably  causing  fracture, 
can  be  calculated  by  marking  off  the  crack  size  on  the  negative- A'  axis  of  an  R-curve,  and 
constructing  the  tangent  to  the  R-curve  passing  through  this  point.  The  slope  of  this  curve 
is  proportional  to  the  load/  fracture  stress  at  which  this  crack  propagates  unstably. 
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crack  size  (microns) 

Figure  6-10:  This  shows  the  fracture  stress,  as  a  percentage  loss  from  the  strength  of 
structures  with  no  cracks,  measured  for  2-,  3-  ,  and  4-hierarchy  structures  for  different  crack 
sizes.  The  shaded  region  shows  the  crack  sizes  for  which  there  is  a  10%  loss  in  strength. 
VVc  observe  that  the  sensitivity  of  fracture  strength  to  crack  size  goes  down  with  increasing 
hierarchy  level. 


structures  from  sub-micron  length  scales  to  several  hundreds  of  microns  in  size.  The 
basic  constituent  of  these  structures  is  silica,  which  are  nanostructured  into  porous 
form,  inspired  by  diatom  morphology,  as  shown  in  Chapter  4.  Composites  are  created 
out  of  nanoporous  silica  and  bulk  silica  as  constituents. 

We  begin  the  study  with  the  design  of  periodic  ordered  composite  structures 
of  bulk  silica/ nanoporous  silica,  classifying  them  into  two  representative  systems; 
one  in  which  the  soft/tough  material  is  the  continuous  phase,  with  the  hard/brittle 
material  dispersed  in  platelet  form  within  the  matrix,  mimicking  several  hierarchy 
levels  in  bone  and  nacre.  The  second  representative  system  is  where  the  hard/brittle 
constituent  is  the  continuous  phase,  with  regions  of  soft/tough  material  embedded 
within  the  hard  matrix.  These  structures  are  found  in  biocrystals,  where  protein 
material  is  encapsulated  in  a  hard  crystal,  such  as  biological  calcite  single  crystals. 
We  observe  that  the  biocalcite-like  structures  have  a  larger  defect  tolerance,  i.e..  the 
sensitivity  in  change  in  fracture  strength  with  crack  size  is  smaller.  We  observe  that 
this  can  be  correlated  to  diffuse  microcracking  throughout  the  material  even  in  the 
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Figure  6-11:  A  defect-tolerant  length  scale  is  calculated  to  be  the  crack  size  at  which 
the  material  retains  90%  of  its  fracture  strength.  This  figure  shows  how  this  length  scale 
increases  non-linearly  with  the  number  of  hierarchies.  The  defect-tolerance  reaches  up  to 
~60  jxm  with  4  levels  of  hierarchy.  The  red  line  represents  an  exponential  fit  of  L  = 
1.17e12D5-v,  where  L  is  the  defect-tolerant  length  scale  in  microns,  and  N  is  the  number  of 
hierarchy  levels.  The  limiting  value  of  1  level  of  hierarchy,  i.e .,  bulk  silica,  does  not  show 
any  defect-tolerant  response  to  crack  size. 


presence  of  stress-concentrators  such  as  small  cracks. 

Subsequently,  we  design  3-  and  4-hierarchy  level  structures  using  different  com¬ 
binations  of  the  bone-like  and  biocalcite-like  templates.  We  observe  that  the  defect 
tolerance  length  scale  increases  with  number  of  hierarchy  length  scales.  R-curves. 
which  capture  change  in  fracture  toughness  per  unit  length  stable  crack  advance  in 
a  material,  arc  measured  for  all  hierarchy  levels.  The  defect- tolerance  property  is 
then  correlated  with  a  rising  R-curvc  behavior,  in  particular  the  slope  of  the  R-curvc. 
It  is  seen  that  the  R-curve  increases  in  both  magnitude  and  slope  as  the  number  of 
hierarchy  levels  increase.  Experimentally,  R-curves  have  been  measured  for  fracture 
samples  of  hierarchical  biological  materials  such  as  bone,  nacre,  and  dentin,  and  a 
similarly  rising  R-curve  behavior  is  seen  [190.  184.  197,  198]. 

The  improvement  in  the  defect-tolerance  property,  which  captures  both  increase 
in  crack  initiation  toughness  [30],  and  the  increase  in  the  R-curve  slope,  is  captured 
numerically  through  a  length  scale  at  which  t  he  structure  loses  10%  of  its  no-crack 
fracture  strength.  We  observe  that  this  quantity  increases  rapidly  with  number  of 
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hierarchy  levels,  reaching  «60  pm  for  4  levels  of  hierarchy. 

The  use  of  several  levels  of  hierarchy  is  thus  shown  to  impact  toughness  through 
many  effects.  Apart  from  higher  absolute  toughness  values,  we  find  that  the  use  of 
hierarchies  also  increases  the  rate  of  increase  of  toughness  per  unit  crack  advance. 
These  together  lead  to  a  higher  defect  tolerance  for  structures  with  greater  number 
of  hierarchy  levels.  This  behavior  can  provide  an  explanation  for  the  use  of  multiple 
hierarchies  in  biological  materials.  Stable  crack  advance  is  not  catastrophic  or  dis¬ 
abling  for  biological  materials,  since  many  materials,  such  as  bone,  have  the  ability 
to  self-heal  over  time.  A  rising  R-curve  behavior,  across  several  microns  or  higher 
length  scales,  promotes  stable  crack  advance  in  the  material  and  allows  it  to  be  useful 
beyond  its  fracture  initiation  load  point.  Thus  large  loads,  which  would  shatter  a 
single  hierarchy  material  with  a  single  crack  propagating  right  through  the  material, 
dissipate  energy  in  multiple-hierarchy  materials  through  several  cracks  being  initiated 
and  arrested  at  different  length  scales  (sec  Figures  6-3(b)  and  6-7(b)  and  (e)).  If  the 
load  doesn’t  rise  to  unstable  propagation  values,  this  probably  provides  the  material 
time  to  heal,  and  still  carry  load.  We  observe  that  a  higher  number  of  hierarchies  is 
provides  a  better  rising  R-curve  behavior. 

Future  work  could  be  aimed  at  optimizing  the  entire  R-curve  behavior  over  several 
microns  or  higher  length-scale  of  crack  advance.  Both  the  number  of  hierarchies, 
and  design  of  individual  hierarchy  levels,  would  have  different  effects  on  the  entire 
R-curve  shape.  Our  mesoscale  model  is  a  stepping  stone  that  can  be  utilized  for 
design-optimization  to  maximize  R-curve  magnitude  and  slope  improvements.  The 
model  can  also  be  extended  to  larger  length  scales,  potentially  sub-nun,  such  that  R- 
curvo  toughness  improvement  is  seen  over  length  scales  comparable  to  the  macroscopic 
sample  size. 
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Chapter  7 
Conclusion 


This  chapter  summarizes  key  findings  of  this  thesis,  and  discusses  the  significance  of 
results  obtained  from  simulation  and  theoretical  considerations  of  the  importance  of 
(a)  nanostructuring  and  (b),  the  use  of  hierarchical  assembly  in  bio-inspired  materials. 
The  ability  to  use  mechanically  'poor5  materials  as  building  blocks  to  obtain  structures 
that  optimize  disparate  mechanical  properties  simultaneously,  is  the  goal  of  mimicking 
biological  mechanical  design,  and  a  bottom-up  design  methodology  to  achieve  this 
has  also  been  outlined.  Also,  later  in  this  chapter,  we  touch  on  directions  for  future 
research  to  improve  and  build  on  these  results. 


7.1  Summary  of  key  findings  and  significance 

Mimicking  biological  materials  design  is  an  important  frontier  in  materials  design. 
This  is  driven  by  several  enviable  properties  of  biological  materials,  primary  of  which 
is  they  utilize  very  ’poor’  engineering  materials  such  as  protein,  silica,  hydroxyapatite 
to  build  animal  skeletons  and  protective  armors.  They  also  combine  these  ’poor’  en¬ 
gineering  materials  to  make  composites  that  possess  a  combination  of  the  best  of 
the  properties  of  their  constituents.  As  an  example,  this  behavior  is  easily  captured 
through  Ashby  materials  properties’  charts  of  stiffness  versus  toughness  for  biological 
materials  and  engineering  materials  (see  Figure  1-1).  Engineering  materials,  such  as 
metals,  alloys,  ceramics,  and  their  composites,  show  a  ‘banana-curve’  behavior  where 
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one  of  these  properties  is  improved  at  the  cost  of  another  [4].  Biological  materials, 
made  of  protein  and  calcit.e /hydroxyapatite ,  such  as  bone,  antler,  enamel  etc.  show 
an  ‘inverse  banana-curve’  behavior  where  these  materials  combine  the  high  toughness 
of  the  protein  constituent  with  the  high  stiffness  of  the  mineral  constituent  [2] .  Bio¬ 
logical  materials  also  present  multi-functionality  and  optimization  of  several  desirable 
properties.  Finally,  they  also  show  robustness,  or  resistance  to  catastrophic  damage. 
Due  to  the  heterogeneous,  hierarchical  nature  of  the  design  structure  of  these  biolog¬ 
ical  materials,  and  the  presence  of  a  wide  variety  of  designs,  it  is  key  to  be  able  to 
discover  key  underlying  universal  design  principles  that  can  then  be  transferred  and 
applied  to  any  common  engineering  constituents.  In  this  work,  we  identified  three 
general  design  principles, 

•  nanostructuring  according  to  nano-motifs  found  in  several  biological  systems, 

•  use  of  ‘poor’  materials  as  structural  building  blocks, 

•  use  of  several  hierarchy  levels  to  reach  from  the  nanoscale  structure  to  macroscale, 
using  well-defined  assembly  at  each  level  of  intermediate  hierarchy. 

We  began  our  analysis  by  firstly,  selecting  the  platelet-matrix  nanostructure  of 
bone  and  nacre,  as  a  design  template  to  mimic.  Since  these  nanostructures  are  com¬ 
posed  of  hard/stiff  calcite/hydroxvapatite  platelets  embedded  in  a  soft/ductile  pro¬ 
tein  matrix,  we  chose  hard  and  soft  metals;  to  replace  the  bone  constituents.  We  thus 
designed  metal-matrix  composites  with  geometric  design  inspired  by  bone  nanostruc¬ 
ture  in  Chapter  3,  and  studied  their  mechanical  properties  under  tensile  loading.  We 
found  significant  dependence  of  flow  strength  of  these  nanocomposites  on  the  geome¬ 
try  parameters  of  the  bone  nanostructure,  and  found  size  limits  which  maximized  the 
strength  of  these  composites.  Or  analysis  also  revealed  that  there  exist  fundamental, 
intrinsic  length  scales  that  control  their  plastic  deformation  mechanism  and  strength 
properties.  Specifically,  we  found  that  the  use  of  elongated  platelets  in  the  composite, 
of  high  aspect  ratio  and  staggered  arrangement  of  platelets  for  optimal  load  transfer, 
and  control  of  spacing  between  layers  of  platelets  are  critical  factors  in  strengthening 
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the  material.  The  spacing  between  platelets  perpendicular  to  the  direction  of  loading 
was  found  to  have  a  Hall-Fetch  like  effect  on  the  flow  strength.  We  also  observed 
that  the  strength  of  the  platelet /matrix  interface  determines  the  optimal  size  of  the 
platelets  at  which  the  maximal  flow  strength  is  observed  for  the  nanocomposite.  The 
interfacial  strength  between  matrix  and  platelets  turned  out  to  be  key  in  determin¬ 
ing  the  importance  of  sliding  and  decohesion  as  deformation  mechanisms  at  small 
platelet  sizes,  and  inter-platelet  spacing.  This  leads  to  a  peak  flow  strength  for  the 
nanocomposite  as  a  function  of  platelet  size.  Analysis  of  the  mechanisms  behind  this 
plateauing  of  strength,  revealed  that  there  exists  a  few  fundamental  length  scales 
that  depend  only  on  material  parameters  and  the  particular  geometry,  that  control 
the  plastic  deformation  mechanism  in  small  crystals  under  confined  conditions.  The¬ 
oretical  analysis  shows  that  these  characteristic  length  scales  separate  regimes  of  no 
dislocation  activity,  partial  dislocation  plasticity,  and  complete  dislocation  plasticity 
at  crack  tips  in  ductile  metals.  Wc,  confirmed  this  effect  by  direct  atomistic  simulation 
of  shear  loading  of  a  model  ductile  single  crystal  system.  Thus,  geometrically  confined 
ductile  phases,  in  the  metal-matrix  nanocomposites  under  loading,  will  show  a  tran¬ 
sition  to  homogeneous  shear  based  plasticity  below  a  critical  length  scale.  This  could 
provide  important  guidance  for  the  optimal  design  of  such  nanocomposites,  as  below 
a  certain  size  scale  of  the  ductile  constituent,  it  will  fail  at  its  theoretical  strength  and 
any  further  reduction  in  the  critical  dimension  will  not  increase  the  failure  strength. 

Metals,  however  are  economically  costly  constituents  to  use  in  design.  Is  it  possible 
to  design  bio-inspired  structures  using  cheap  and  readily  available  materials?  An  ex¬ 
ample  is  silica,  which  is  abundantly  found  as  sand:  economically  very  inexpensive  but 
also  possessing  ‘poor1  mechanical  properties.  In  particular,  it  is  brittle  and  possesses 
low  fracture  toughness,  and  cannot  be  used  as  a  structural  material  in  its  bulk  form, 
whether  its  amorphous  glassy  form  or  crystalline  quartz  form.  However,  in  biological 
materials,  silica  is  abundantly  used  to  create  protective  exoskeletons  in  many  species 
of  diatoms  and  sea  sponges.  We,  thus,  next  moved  to  the  use  of  this  ‘poor’  engineering 
material,  silica,  instead  of  metals  for  creating  bio-inspired  structures  in  Chapter  4. 
We  turn  to  the  structure  of  diatom  exoskeleton  (see  Figure  1-2)  for  inspiration.  The 


151 


mechanics  of  diatom-inspired  nanoporous  silica  structures  under  tensile  loading  were 
investigated  using  molecular  dynamics  simulation  and  theoretical  analysis.  We  found 
that  the  elastic  modulus  of  the  nanoporous  structures  is  geometry  dependent,  and  can 
be  modified  for  a  given  porosity  by  changing  the  microstrueturc.  Moreover,  we  found 
that  the  stress-strain  loading  curves  for  these  silica  structures  exhibit  plasticity  be¬ 
low  a  certain  size  scale.  We  were  able  to  establish  a  size-dependent  brittlc-to-ductile 
transition  in  nanoporous  silica  structures.  The  aspect  ratio  and  shape  of  pores  can  be 
modified  to  control  the  yield  strength  of  these  structures,  and  obtain  large  amounts 
of  ductility  of  up  to  120%  strain.  The  structures  that  show  plastic  yielding  also  show 
large  toughness  improvement  over  bulk  silica.  These  results  reveal  that  nanostruc¬ 
turing  with  control  of  porous  geometry  can  lead  to  application  of  silica  in  carrying 
loads  in  small  devices. 

The  advantages  of  this  nanoporous  design,  however,  come  at  a  cost;  these  struc¬ 
tures  lose  up  to  90%  of  the  stiffness  of  bulk  silica.  Since  hierarchical  design  is  hypoth¬ 
esized  to  be  a  methodology  to  achieve  multiple  property  optimization,  is  it  possible 
to  use  hierarchies  to  improve  stiffness  of  the  nanoporous  silica  structures  up  to  bulk 
silica  values?  Moreover,  can  this  be  accomplished  by  engineering  silica  alone,  i.e.,  use 
of  a  single  constituent  material?  In  Chapter  5,  we  demonstrated  that  it  is  possible  to 
enhance  the  toughness  of  silica  while  retaining  its  stiffness,  without  using  any  other 
material,  through  the  use  of  structural  hierarchies.  We  developed  an  atomistically- 
informed  mesoscale  model  that  can  model  bulk  silica  and  nanoporous  silica  behavior 
at  the  micron  length  scale.  By  design  of  a  randomly-dispersed  composite  structure 
of  the  bulk  silica/  nanoporous  silica,  we  are  able  to  retain  stiffness  up  to  70%  of 
bulk  silica  while  increasing  toughness  four  times  over  its  value.  This  improvement 
in  toughness  and  retaining  of  stiffness  of  bulk  silica  by  designing  small  regions  of 
nanoporous  geometry  (f»14%)  in  the  bulk  phase  point  to  a  design  methodology  for 
obtaining  stiff  and  tough  materials  out  of  an  inherently  brittle  material  (silica)  to 
begin  with.  Thus  we  were  able  to  achieve  a  hierarchical  structure  with  high  stiffness 
and  toughness  through  the  use  of  a  single  material  silica  which  is  traditionally 
considered  an  undesirable  (mechanically  inferior)  structural  material  due  to  its  brit- 
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tleness,  by  arranging  it.  in  a  hierarchical  pattern  with  substructures  that  scale  from 
the  nanoscale  to  microscale  dimensions. 

We  then  proceeded  to  the  analysis  of  several  levels  of  hierarchical  arrangement  of 
bulk  silica/  nanoporous  silica  composites  on  toughness,  in  Chapter  7.  We  observed 
that  the  structural  11-  curve  shows  a  rising  characteristic  for  additional  levels  of  hi¬ 
erarchy.  Thus  both  the  absolute  value  of  toughness,  and  the  slope  of  the  R-curve 
increase  with  the  number  of  hierarchy  levels  (see  Figure  6-8.  This  leads  to  an  increas¬ 
ing  defect-tolerance,  i.e.,  fracture  strength  loses  its  sensitivity  to  changes  in  crack 
length  (sec  Figure  6-10).  We  showed  a  non-linear  increase  in  this  defect-tolerance 
length  with  additional  levels  of  hierarchy  (sec  Figure  6-11).  This  behavior  can  be 
linked  to  the  use  of  multiple  hierarchies  in  biological  materials.  Stable  crack  advance 
is  not  catastrophic  or  disabling  for  biological  materials,  since  many  materials,  such  as 
bone,  have  the  ability  to  self-heal  over  time.  A  rising  R-curve  behavior,  across  several 
microns  or  higher  length  scales,  promotes  stable  crack  advance  in  the  material,  and 
allows  it  to  be  useful  beyond  its  fracture  initiation  load  point. 

This  thesis  established  a  direction  and  methodology  in  understanding  and  applying 
universal  design  principles  that  can  be  gleaned  from  biological  materials  for  materials 
design  for  structural  applications.  Fully-atomistic  and  mesoseale  modeling  showcased 
a  bottom-up  fully  computational  approach  to  implementing  these  design  principles. 


7.2  Opportunities  for  future  research 

The  previous  section  discussed  the  atomistic,  theoretical  and  mesoseale  modeling  of 
implementation  of  bio-inspired  design  principles  in  hierarchical  structures  based  on 
metals  and  silica.  The  purpose  of  this  section  is  to  illustrate  the  shortcomings  of  the 
current  work  and  propose  directions  for  future  research  on  computational  design  of 
bio-inspired  hierarchical  materials. 

The  design  and  manufacturing  of  bone-inspired  metal- matrix  nanocomposites  re¬ 
quires  the  simulation  of  actual  metal/metal  composite  systems.  Model  materials,  as 
used  in  Chapter  3,  provide  broad  design  guidelines,  but  the  question  of  exactly  which 
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metal  constituents  to  use  can  be  settled  by  atomistic  simulations  with  well-developed 
metal  and  metal  alloy  interatomic  potentials.  The  thermodynamic  and  kinetic  sta¬ 
bility  of  interfaces  between  the  constituents  and  the  equilibrium  configuration  of  the 
interfaces  are  also  important  considerations  to  the  mechanical  behavior  that  can  only 
be  settled  through  simulations  of  particular  metal/metal  systems  and  experiments. 

The  study  of  the  design  of  diatom-inspired  nanoporous  silica  structures  can  be 
enhanced  by  analyzing  their  mechanical  properties  in  a  hydrated  environment.  This 
is  because  the  presence  of  water  has  been  shown  to  have  an  effect  on  silica  mechanical 
properties,  both  in  bulk  and  nanoscale  [157].  Some  preliminary  investigations  have 
been  performed  on  the  response  of  these  nanoporous  structures  in  the  presence  of  hy¬ 
drogen  and  water  [156].  Silica  surface  termination  is  also  affected  by  the  environment 
and  nature  of  the  exposed  crystal  surface,  and  the  nanoporous  structures  possess 
sufficient  surface  area  such  that  the  surface  mechanical  response  is  important.  Also, 
these  designs  experiments  can  be  undertaken  for  amorphous  silica,  which  is  the  form 
in  which  bio-silica  in  diatoms  is  found. 

The  study  of  hierarchical  assemblies  of  silica  through  atomistically-informed  mesoscale 
models,  as  seen  in  Chapters  5  and  6,  can  be  regarded  as  a  preliminary  attempt  in  this 
emerging  area  of  computational  study  of  hierarchies  in  materials.  Mesoscale  models 
have  been  used  for  several  other  biological  systems  with  hierarchies,  such  as  collagen 
fibrils,  amyloids,  spider  silk,  and  is  a  promising  approach  to  extracting  qualitative  de¬ 
sign  information  about  the  requirement  and  advantages  of  hierarchical  design.  This 
can  be  contrasted  with  fractal  continuum-level  approaches  [36,  187]  which  arc  much 
coarser  in  the  details  of  geometry  and  deformation  they  capture.  In  further  studies, 
mesoscale  models  of  hierarchical  materials,  can  be  used  along  with  design  optimiza¬ 
tion  techniques  such  as  genetic  algorithms,  to  optimize  material  distribution  and 
arrangement  over  several  length  scales  of  hierarchy.  In  particular,  it  may  be  possible 
to  optimize  R-curve  shape  over  several  levels  of  hierarchy  and  length-scale  to  increase 
crack  propagation  toughness,  and  stop  and  repair  cracks  before  they  cause  complete 
failure.  \ 
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